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That mathematics plays a certain role in various phases of economic theory is, of 
course, quite well known. The number of mathematicians and economists unaware of 
the discipline known as ‘““Mathematical Economics” is surely small. Unfortuntely, the 
number of mathematicians and economists who are aware of the aims, the methods and 
the successes (the failures are, as is always the case, only too well known) of this rapidly 
expanding field, is likewise small. 

The majority of us are exposed in the very early stages of our training to simple 
applications of mathematics in economics. To cite an instance, the classical high-school 
course in algebra which does not use the price-demand relation to illustrate the ideas of 
an inverse variation, is probably non-existent. A good many mathematicians think that 
this is almost the sum total of pure mathematical depth and sophistication that the 
economist encounters. 

On the other hand, since statistics is, in so many obvious ways, ideally suited for 
analyzing certain types of economic phenomena, one is often fooled into believing that 
this represents all of the applied mathematical ideas that can play a central position in 
an economic investigation. 

The appearance, in recent years, of the book “The Theory of Games and Economic 
Behavior”? by von Neumann and Morgenstern [40] has been instrumental in dispelling 
from the minds of many mathematicians and economists some of the false ideas about 
what mathematics entered, and possibly more important, what mathematics does not 
enter, in the problems arising in everyday mathematical economics. 

The purpose of this paper is, in part, to give a presentation of some phases of pure 
mathematics that are in current use in the economic world. We do not claim or pretend 
that this paper is either exhaustive or definitive—we merely propose to touch on a few 
boundary points with the hope that some of the readers will feel an urge to dig deeper 
into the interior. This paper is being written by a mathematician primarily for a math- 
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ematical reader. The topics we discuss are thus selected, not so much for their economic 
depth or significance as for their mathematical interest. 

Historically, the calculus entered very quickly into the study of price-demand re- 
lations. The calculus was, in fact, essentially the only tool used, until fairly recently in 
all of mathematical economics. (For a history and description of both the methodology 
and problems considered up to the beginning of the 20th Century see the authoritative 
article by V. Pareto [43]). 

To say that the calculus, as a fundamental weapon, has by now played out its role 
would be both exaggerated and misleading. However, its initial foothold has been weak- 
ened considerably, and in certain places, arguments previously employing calculus 
techniques have been reworked using more powerful and more modern methods. There 
are several fundamental reasons for this trend to get away from using the calculus. 
Firstly, in a desire for generality, the conditions of differentiability placed on functions, 
especially when not needed, are aesthetically unsatisfactory. Secondly, the existence of 
derivatives for functions arising from a study of complex and oftentimes highly discrete 
situations is by no means a natural assumption. Thirdly the appropriate calculus con- 
ditions imposed on the functions considered at times obscure the essential nature of the 
problem. Added to this are the huge successes achieved in other fields (physics, for one) 
by the introduction of more of the full mathematical power available today. 

This accounts for the tendency, with certain economists, to introduce some ‘‘natural’”’ 
tools in their domain. The surprising effect (to some people) has been a great simplifica- 
tion and reinterpretation of old results and a satisfying surge forward in new research. 

Although it is clear that mathematical statistics and game theory play a vital part in 
economic theory, we shall not consider their applications here. How and why they are 
used can be found very easily, and in many places [3, 33, 38, 40, 41, 50). 

As we have pointed out previously, the calculus has performed (and is still per- 
forming) a striking function in the discipline. We begin the paper proper with an example 
that arises in economics and is handled via the calculus. This example has its origin in 
the theory of consumer’s behavior; the derivation we present is of the so-called Slutsky 
Equation (Slutsky [47]; for a recent treatment of it see Samuelson [45, p. 97—103]). 
The problem itself is of classical stature in the economics of consumption. Its solution is 
a striking example of how a serious economic result emerges from an elementary math- 
ematical analysis.” 

We suppose there are n commodities which are labeled 1, 2, --- , m and a given 
consumer. A commodity bundle is a vector x whose ith coordinate represents an amount 
of the ith commodity. We suppose that we are also given prices for each commodity. 
Let p be the price vector. Given a commodity bundle z it is assumed that the consumer 
derives a certain satisfaction s(x) where s is a real-valued, twice differentiable strictly 
convex function. (Economists usually consider a more general class of function, but for 
the sake of simplicity we restrict ourselves to the convex case). For economic purposes s 
could be equally well replaced by a monotonic increasing transformation of itself. In 
order for our results to be economically meaningful, our description of the consumer’s 
behavior should be independent of the particular transformation applied to s. Finally, 
we suppose the consumer has a given amount of money, u. Subject to the budgetary 


*The form of the proof given here is due to G. Debreu. 
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constraint 
px =u (’ denotes transposes) (1) 


the consumer tries to maximize his satisfaction s(x). Out of this desire, there results 
the necessary condition 
ds 
dx 


where oc, of course, depends on the particular form of s. (¢ has traditionally been called 


= op, o a scalar and (2) 


the marginal utility of money). 
On differentiating (2) we obtain 


Sdx = adp+pde (3) 


where S is the Hessian matrix of s. S is thus a symmetric matrix. 
Differentiation of (1) yields 





p’ dx + (dp')x = du. (4) 
Let 
ts r| 
7 - |” 
p’ 0) 


From the symmetry of S, we have that = is also symmetric. The system (3) and (4) 


now becomes 


dx ( dp | 
z = | (5) 
- l ay —z’ dp) 


40 ° . —=1 . 4° — s . 
Since s is strictly convex, =~ exists. Since = is symmetric, 
[x 7 
, 

/ 
—— 


where X is again symmetric and where y is a column vector and c is a number. Equation 


gt: _. 


_ — 


(5) then gives rise to 


| dx | (x . ( dp | ‘i 
—_ , 6 
|- Ly” c m —x’ dp) 


Expanding the right-hand side of (6), we arrive at dr = X dp + vy (du — x’ dp). Thus 


eae Se 
Op vs ad Ou 
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Combining these we have 


Ox 2 Ox ., ss 
a Oe ea 2 2 (4) 
op Ou 
and so 
: Ox OL, 
=> .—— 2. (3) 
Op Ou 


Since X is symmetric, picking out the (7, 7) and (j, 7) element we have 
Ox Ox Ox Ox 
—— 2 a oe 2 (9) 
Op; ' Ou Op Ou 
This is the Slutsky Equation. Notice that the result is independent of the particular 
transformation applied to s. The respective sides of (9) are termed the substitution co- 
efficients. We now seek an economic interpretation for these. We ask: what changes in 
amounts of commodities and what changes in income would leave the consumer’s satis- 


faction unchanged. From (2), since s is being kept constant, p’ dr = 0. So (4) becomes 
du = x’ dp. Thus dx = X dp, whence the matrix elements of X, that is, the substitution 


coefficients describe the consumer behavior when satisfaction is assumed to be held 
constant. 

This simple illustration of the use of the calculus, and of a calculus of the most rudi- 
mentary form, accounts, to a small degree, for the security the economist felt in the con- 
tinued use of the calculus as his primary tool. As a consequence, he felt no particular 
urgency to broaden the mathematical repertoire to be applied to his problems. I should 
not like to imply that the preceding example is typical, in its mathematical depth, of the 
category of problems attacked. Mathematical analysis of a far more penetrating character 
was employed by some economic writers. Nontrivial existence theorems in the theory of 
differential equations, results from the theory of integral equations, and a great wealth 
of related mathematics was put to very effective use. (For an excellent treatment, along 
such lines, of a wide variety of topics, see Samuelson’s stimulating book [45].) 

However, in the late 1930’s, the variety of mathematical topics finding application 
became more and more diversified and modern mathematical developments began re- 
flecting on the nature of the work undertaken by the economist. In the remainder of this 
paper we shall touch on some of these. The order of their presentation is not meant to be 
chronological. These often involve convex sets or the calculus; this is quite’ natural, 
since such an important part of economics can be viewed as a maximizing activity. 

We now consider a situation arising in the phase of economics known as “Welfare 
Economics.’”’ Before discussing the problem per se we need a brief, crude description of 
the general scope of the Welfare Economics. (For a much more thorough description see 
Samuelson [45], Arrow [1]). 

Suppose that there are n-individuals in a community, designated by 1, 2, --- , n. 
We suppose that each of these has a preference ordering by which he ranks the possible 
prospects or social states that can confront him; that is, to the ith individual there 
belongs a binary relation R; defined on the set of possible social states X, Y, Z with the 
following properties: 


1) XR,Y or YRX, 


2) XR, Y, YR;Z implies XR;Z. 
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Here, R; is a complete (or as it is sometimes called, a weak simple) ordering; XR;Y can 
be thought of as ‘‘X is at least as good as Y as far as the ith individual is concerned.’”’ 
Let XP,Y denote that XR,;Y but not YR;X. A social state X will be said to be optimal 
if the following is true: if for any 7 there exists a Y with YP,;X, then for some other 
individual 7, XP;Y. That is, X is optimal if whenever a social state Y is preferred to X 
by some individual, then Y is not ‘at least as good as X” according to all other indi- 
viduals. Another way of looking at this type of optimality can be formulated as: let 
XRY mean XR,Y for all i = 1, 2, «++ , n; R then defines a partial order in the set of 
social states; a state X is now optimal if it is a maximal element in this partial order. The 
problem of the Welfare Economics is to give a description of (and a prescription sa to 
how to attain) such an optimal social state. 

Here is a situation which has been handled within the last ten years by two totally 
different approaches—one in the framework of the orthodox calculus and the other in 
that of the theory of convex sets. In order to exhibit the differences (and the metamor- 
phosis) in these two polar attacks on the problem we present two solutions of the same 
problem arising in the Welfare Economics. The first treatment we give is the calculus 
discussion, and is due to Lange [35]. 

We assume each individual’s preference ordering R; can be represented by a real 
valued function u“”, the ith person’s utility function. The purpose is to maximize u“’’, 
for each i, subject to the conditions that u‘” is kept constant for 7 ¥ 7. 

We suppose the utility function of each individual to be a function of the amounts of 
each commodity that he gets; that is, if x"? is the amount of the jth commodity held 
by the 7th individual, 7 = 1, 2, --- , m, then u’? = u“?(x;"’, 2," --- , x’). We further 
assume that each u‘” is appropriately differentiable. Let X, = >.7-, 2‘"’ be the total 
amount of the rth commodity which is obtained by the whole community. The X, are 
interrelated by a “technological function” F(X, , X., --: , Xn) = 0. The problem then 
becomes: maximize u‘‘’(x;"’, --- , x) t = 1, 2, «++ , n subject to the constraints 


1) u’?’(ay"’, +++ , 2m’) = constant for j # 7, 


I 


2) X, = Do 2x!°" for each r = 1, 2, --+ , m, 
i=1 
3) F(X, , X2,°°: , Xn) = 0. 
From the theory of Lagrange multipliers, this is equivalent to “‘extremizing”’ the 
unconstrained expression 
1) r 7 r . 
Yaa? + Sn Sa? — X,) + K+ Xd, 
l r=1 i=1 
where the \’s and yw’s are Lagrange multipliers. Differentiation and elimination yields 


a" Fae” oF OF . 
. Sa | for all 7, r and s. 
Oz,’ ox,"” on. @2. 


Equivalently, 


au? Jadu®? — au? Sau“ . - 
— = ——— ' for all 7, j, r and s. 


a .(¢) ( ee ( o 
ox*" ore” asf axis” 
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Economically, du‘’’/dz,"’ is the “marginal utility of the rth commodity for the ith in- 
dividual,” and the above gives, as a necessary condition, that certain ratios of marginal 
utilities be equal. The result is invariant with regard to monotone increasing trans- 
formations of the utility functions. The Lagrange multipliers could be interpreted as 
prices, and the solution of the problem could be stated in terms of the existence of prices 
with certain properties. This will be done in the second treatment we give of this same 
problem using the theory of convex sets. 

Before turning to the second version of this problem, we digress into several path- 
ways suggested by factors which have arisen in the examples which we have already 
discussed. 

In the situations which we have considered, the following kind of assumptions have 
been made: a complete ordering, occurring naturally as a preference ranking, confronts 
the economist; he, in turn, wishes to discuss the optimal behavior under this ordering. 
In order to do so he resorts to a real-valued function, having desirable differentiability 
properties. However, the final description of the optimal behavior is given in a form 
independent, to a major degree, of this function. Two questions immediately present 
themselves. Firstly, what conditions on the ordering insure the existence of such order- 
preserving functions? Secondly, since the final form of the solutions does not depend on 
these functions, can all these problems be handled without the artifice of a real-valued 
representation of the ordering? 

The first of these problems, important as it seems to be, has received very little 
attention from the economist. In fact, some economists have even gone as far as to 
tactily assume that every complete ordering can be so represented. (The lexicographic 
ordering of the unit square offers an easy counter-example to this). In case the completely 
ordered set is, in addition, endowed with a “probability-mixing”’ operation, this problem 
has been thoroughly thrashed out by von Neumann-Morgenstern, Marschak [37], and 
Herstein-Milnor [25]. For the case in which no such mixture operation exists, H. Wold 
[54] did give conditions on the ordering which guaranteed the existence of an order- 
preserving function; however, his conditions were somewhat restrictive. In the math- 
ematical literature, in a paper by Eilenberg [15] there is given a fairly general set of 
conditions for the existence of a continuous order-preserving function. This has recently 
been rederived by Debreu [10]. The result can be stated as follows: 


Let S be a completely ordered topological space such that 


I) S is separable and connected 

II) for every a) e S the sets fae S| a S ap}, {ae S| a S a} are closed. Then there 
exists on S a continuous, real, order-preserving function. 

Since the economist usually works in a finite-dimensional Euclidian space, or at worst, 
in a separable Hilbert space, Debreu’s condition that his space be separable and con- 
nected is really not very restrictive. Incidentally, Debreu has recently shown that if 
(I) is replaced by: ()’ S is perfectly separable, then the final result still holds. 

As for the second problem, namely that of entirely avoiding the real-valued repre- 
sentation of the ordering, the alternative approach, which we are about to present, to 
the problem in Welfare Economics furnishes a typical method which has been success- 
fully employed. Although this was hardly the main goal of the paper to be cited, it was 
an interesting by-product of the approach. The discussion which follows stems from a 
paper by Debreu [11]. 
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Let us again suppose that there are m commodities labeled 1, 2, --- , m and n con- 
sumption units, labeled 1, 2, --- , n. The activity of a consumption unit is characterized 
by a vector x, the consumption vector, in the m-dimensional Euclidian space; the ith 
component of this vector being the amount of ith commodity used by the consumption 
unit. We suppose that each consumer has his individual complete order R; defined on 
the set of all consumption vectors. It may happen that zR,y and yR,x without y = 2; 
but if we define x/,y if this does occur, this yields a proper equivalence relation, and the 
equivalence classes, S,;(x), now form a linearly ordered set 7; . The satisfaction space $ 
is introduced as the set of all vectors S = (S,, S,,---, S,) where S;e7; . In § a partial 
order is defined by 


¥(2) (2 ¥(2) (2) (1) (1) cy(1) (1) 
S = (8; , Se -- S)2S8 = (S; » os, °°, &) 


1 
if and only if S;”R;S;" for every i = 1, 2, ---,n. 

The production activity of the system is represented by an input vector y = (y% , 
Yo, °** » Ym), Where y; is the input (positive or negative) of the 7th commodity. Tech- 
nological conditions restrain y to belong to a set 7. Let n”’” be the set of efficient produc- 
tion vectors. A family of sets n; ,j7 = 1,2, --- , risa decomposition of nif 7 = } » n; , the 
sum being in the sense of vector sums of sets. 7; characterizes activities of the jth pro- 
duction unit. If X¥; = X,(s?) denote the {z, | s;(z;) = s'} the ith consumption unit, 
X = >> X, is the set of all total consumption vectors. z = x + y, x eX, y « 7 is the total 
net consumption of the whole economy. Let z° be the vector whose components are the 
available amounts of each commodity. The economic system is constrained by 


0 


yen, z2S2 in the vector sense. 


The goal of the economic system is to find an S e S$, maximal in the partial order 
defined on $, which is consistent with the above constraints. We now assume all the 
sets X; , n; are convex and closed. 


Z=>X.+ Da is then convex. 


Using the separation theorem for convex sets, namely, if two closed convex sets with 
interior points have only one point in common then there is at least one plane through 
that point separating the convex sets, one readily obtains the following result: 

A necessary and sufficient condition for S° to be maximal (or for 2. = >, 2? + >> y? 
to be minimal) is the existence of a price vector p > 0 and a set of numbers a; , 7 = 1, 2, 


, n so that 


a) x, being constrained by p’x; = a; , s;(z;) reaches it maximum at zx! for every i 


8) y, being constrained by y; 7; , p’y; reaches its minimum at y for every j. 


Here p is a positive, normal vector to the separating plane. The theorem restates 
the following rules of behavior for consumption and production units: each consumption 
unit, subject to a budgetary constraint maximizes its satisfaction and each production 
unit, subject to the technological constraints, maximizes its profit. 

We leave these problems concerning convexity, maximization under constraints and 
related topics and turn to a phase of economics which employes completely different 


mathematical techniques. 
Matrix theory, to some degree, enters into many portions of many diverse fields. 
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Matrices are used in these as a pure notational device, as a compact and transparent 
representation of systems of linear equations, and in many other subsidiary, convenient 
roles, without use being made of essential matrix theory. In economics, also, matrices 
make their appearance in this costume. But deep, significant, intrinsic results of the 
matrix theory do play an important part. In the recent literature one finds a number of 
research papers in economics which employ certain corners of matrix theory in the 
fullest. (To cite a few examples, we refer the reader to Chipman [9], Metzler [39], Goodwin 
23], Simon and Hawkins [24] and Solow [48].) 

The mathematical origin of a good many of these results is a sequence of papers by 
the German mathematician Frobenius [19, 20]. These results have been rederived and 
somewhat extended in a simple, modern fashion by Wielandt [53] and Debreu and 
Herstein [13]. In these above-mentioned papers, among other results, are obtained 
theorems concerning the nature of the characteristic roots of a non-negative matrix A 
and the properties of (sJ — A)~* for s a sufficiently large real number. 

We now consider an economic situation in which these theorems function effectively. 
We follow here the treatment of a problem in the theory of international trade as given 
by Solow [48]. 

Suppose we consider ” countries carrying the labels 1, 2, --- , n. Let a,; denote the 
marginal propensity (that is, the increase in imports from country 7 by country j per 
unit increase in income of country 7) of the jth country to import from the 7th country, 
and a;; the marginal propensity to consume domestic goods. Furthermore, let x; represent 
the national income of the ith country, c; the autonomous expenditure in country 7. 
Then we have, assuming linearity (which can be viewed as a first approximation for more 


general cases) 
1) Go, = ps a;;«; + ¢, 


In matrix form this becomes 
11) (I -—-az= C; 
where a = (a;;); x, ¢ are column vectors. Economic meaningfulness demands that the 
quantities x; , c; be all nonnegative. We assume that all a;; = 0. 
Thus we are immediately forced to consider conditions under which the systems 
xz = (I — a)" cis solvable in nonnegative terms. 
The linear equation system (ii) is easily seen to be the static solution of the linear 


difference equation system 
ili) Izx(t + 1) — ax(t) =. 


A question of importance for this dynamic system will naturally concern the stability 
properties of its solutions. These can be shown to be equivalent to the condition that all 
the characteristic roots of a are less than 1 in absolute value. For this system Metzler 
[39] has proved: if a;; = O it is necessary and sufficient for the stability of (iii), that 
I — a have all its principal minors positive. (As a pure matrix theorem this is an easy 
consequence of the results of Frobenius, when we use that stability is equivalent to the 
condition that characteristic roots of a be less than 1 in absolute value.) These types 
of questions are treated in Section 3 of the paper by Debreu and Herstein [13]. 
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Let us consider a = (a;;), a nonnegative, n X n matrix. A set, S, of indices will be said 
to be closed if a,, = 0 forge S, pg S. That is, a closed set is associated with a collection 
of countries each of which spends in no country which is not in the collection. If no such 
proper closed set exists the matrix a is said to be indecomposable. This definition co- 
incides with the purely mathematically motivated definition given by Frobenius [20], 
namely: the matrix a is said to be indecomposable if, for no permutation matrix 7, the 
product zar~* can be represented in the form 


(Ain Ais] 
| 
| 


Lo Ass) 


where the A;; are square submatrices. A matrix which is not indecomposable is called 
decomposable. If b is a decomposable matrix then there is a permutation matrix 7 so 
that 


’ 


where the B; are square indecomposable submatrices on the diagonal. 
In terms of the economic model, if the import matrix a is indecomposable, then a 
dollar spent in any one country will eventually induce spending in every other country. 
Let us now assume that the import matrix a is indecomposable. From the results of 
Wielandt [53] or Solow [48] we can find a permutation z which puts a into the form 








0 0 tee 0 G,, 

ee 8 os 0 0 
rar ' = |0 To vee 0 0 |, 

0 0 G.-: 0. 


where the 0’s on the diagonal are square matrices. Economically this can be interpreted 
as: the countries in G; spend only in those in G;,, . 

From the Frobenius theorems the following main result can be extracted: “if a is a 
nonnegative, indecomposable matrix then a has a positive characteristic root r so that 


1) If @ is any other characteristic root of a then | a | S 1; 
2) to r can be associated a positive characteristic vector; 
3) ris a simple root; 

4) an increase of any element of a yields an increase in r. 


Moreover if there should be exactly k roots, a, with | a, | = r then a, = rexp (2r\/k), 
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A = 1,2, --- , k and a permutation 7 exists so that 
0 0 vee 0 ~ 
A, 0 ve 0 0 | 
rar = |0 A» 0 0 | 
(0 0 Se 0 


Solow’s w, which had a purely economic origin now can be interpreted as k, the 
number of characteristic roots of largest absolute value. After w time units have elapsed, 
the original expenditures from G; in G;,, has an influence in G; , giving rise to a cyclicity 
of spending. 

The above discussion encroaches on a tiny part of this territory, which has already 
been staked out by the economist. However, it is typical in the sense that the theorems 
of Frobenius play such a central role. 

We now concern ourselves with several questions which have a combinatorial flavor 
to them. While these are not all phrased purely in economic terms, they have, or should 
have, many applications in economics. 

The first of these is the so-called “personnel assignment”? problem. Suppose that 
there are n individuals available to fill n positions. Suppose further that the 7th individual 
obtains a rating of a;; in his ability to handle the jth job. The question then is: how shall 
individuals be assigned to jobs in order to have the “overall efficiency” a maximum? 

The problem can be put into purely mathematical garb in following vein: given a 
matrix A, for what permutation matrix P is the trace of AP a maximum? Of course this 
would merely require testing n! possible permutations. Even for n relatively small this 
straightforward procedure would be out of the realm of practical feasibility. 

The problem then becomes one of reducing the number of necessary computational 
steps. This can be achieved by (at least) two continuizations of the problem. We shall 
describe one of these, due to von Neumann [41] in very little detail, later in the paper. 
But before doing so, we should like to describe several related combinatorial questions. 

A situation very similar to that in the personnel assignment problem is the ‘‘desk- 
cabinet”’ problem. This runs as follows: Suppose that there are n desks and 7 filing 
cabinets. Let d;; be the distance from the ith desk to the jth cabinet. On the assumption 
that the individuals assigned to each desk make the same number of trips per day to 
their respective filing cabinets, what assignment of desks to cabinets should be made to 
insure that the total distance walked is a minimum. Formulated in mathematical terms, 
the problem is that of finding a permutation matrix P, which minimizes the trace of DP 
where D is the matrix of distances. 

A purely economic variant of this can be phrased as follows: minimize the cost of 
production of n plants at n locations if d;; is the cost of production for the ith plant at 
the jth location. 

A generalization of the two above-mentioned situations, arising fairly naturally in 
economics, involves finding a permutation P, which minimizes or maximizes the trace 
of DP where P is no longer left free to roam over the whole symmetric group but is 
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restricted to range only over a subset thereof. In an actual problem this might be realized, 
for instance, by the condition that in any reshuffling of locations, railway depots would 
be restricted to occupy space only along railroad lines. 

A solution to this problem would immediately yield a solution to another one, the 
“travelling salesman” problem. The travelling salesman problem, verbally described is 
the following: suppose a travelling salesman on his route must visit n cities and return 
to his home base which is one of these cities. If he must visit each of these cities once, 
and if he knows the distances between all of them, how should he plan his trip so that 
the distance he travels is minimal. This can be shown to be equivalent to minimizing 
the trace of CP where P is restricted to be a permutation which in its cycle decompostion 
is representable as an n cycle. 

This relative of the Hamiltonian game can also be exhibited as a special case of a 
purely economic consideration due to Beckmann. Let there be n plants, and n fixed 
locations where these plants could be situated. Suppose a,;; represents the flow from 
the ith plant to the jth plant. Suppose also that k;; is the cost of transportation from 
location 7 to location 7. Then how should the plants be located in order that the total 
transportation costs among plants are a minimum. If A = (a;;) and K = (k;,) and if A’ 
denotes the transpose, this simply becomes a question of finding a permutation P, which 
minimizes the trace of A’P’KP. 

We return to the personnel assignment problem, and present a sketch of von 
Neumann’s game theory approach to it [41]. 

Consider the following two person game: we have an n X n checker board, with 
each square having two indices, its row index and its column index. The first player 
picks a square; the second player then guesses either of the indices of the square which 
the first player has picked. He must state which index he is guessing. If he guesses 
correctly he receives an amount a;; , where 7, j are the indices of the square involved, 
from the first player. Otherwise he receives 0. 

This game is related to the solution of the personnel assignment problem via the 
following theorem proved by von Neumann, (where the strategies refer to those of the 
first player). 

“The extreme optimal strategies of the above game are precisely the following ones: 

Consider those permutations P, which maximize the trace of AP, where A = (1/a;;). 
For each P, assign the probability z,; to the square 7, 7 where x;; = (a/a;;)5p,,;),; where 
6 is the Kronecker delta; and where a is the value of the game for the second player. 

Using techniques of Brown and von Neumann [7] to get approximate solutions of 
games, the number of steps in solving the problem is reduced from n! to a power of n. 

Let z = (z;;) be a vector in the n’ Euclidian space. We define: 


R=(=G@)la20 Deaw=l Lay =, 


f. 
| 
2 

ll 
~ 
x 
nn 
IV 
© 


Be81, Yeas sti, 


P = {z = (@;;) | Ze; = Spey , j for some permutation P}. 


Von Neumann’s proof then hinges on the following two lemmas, 


1) {S = |z z S w forsome we R}, 
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where z < w means the inequality is true in each component; 


2) R = convex huli of P. 


The reader will notice that certain very important current research in economics 
have scarcely been mentioned. For instance nothing has been said about the input- 
output models of Leontief or the applications of linear programming to economics. We 
felt that to do justice to these two would require much more sp::ce than would be appro- 
priate for an article such as this. The reader is referred to Leontief’s book, The Structure 
of the American Economy, |36] for the first topic, the book, Application of Linear Pro- 
gramming to the Theory of the Firm, by Dorfman, [14] for the second topic, and the book, 
Activity Analysis of Production and Allocation [29] for both. 

With this discussion of the combinatorial problems we conclude the phase of the 
paper of “detailed”’ discussion of the applications of mathematics to economics. However, 
before concluding we wish to point out briefly some sources where one can find fine 
applications of other mathematical techniques. 

The theories of differential equations and of difference equations play fundamental 
parts, in many connections, in economic theories. Chapter X of Samuelson’s book [45] 
gives splendid illustrations of the use of techniques from these regions of mathematics. 

Beckmann [5] has employed the classical theory of the calculus of variations in his 
study of continuous models of transportation. Earlier applications of the calculus of 
variation occur in papers by Hotelling [27] and Roos [44] and Evans [17]. 

Stone has used the theory of graphs in his economic studies. An example of this is 
[49]. Fhe theory of graphs has also entered into the considerations of Koopmans and 
Reiter [31] of their transportation model. Charnes [8] has utilized the graph theory to 
extend certain computational techniques. 

Foliowing the trend towards axiomatization in mathematics there have been some 
purely axiomatic studies of economic questions; a pioneer effort in this direction is a 
postulational study of utility by Frisch [18]. Another example is the study of the exist- 
ence of a social welfare function by Arrow [1] and Hildreth [26]. 

Many other applications of the theory of convex sets can be found than those already 
described in this paper. The theory of convex polyhedral cones was first developed by 
Wey] [52], and detailed mathematical investigations of these cones were carried out by 
Gerstenhaber [22]. These results, and many others, have a variety of applications, as 
the reader can see by looking at Koopmans [30], Georgescu-Roegen [21], Samuelson [46], 
and Arrow [2] in the monograph Activity Analysis of Production and Allocation. [29]. 

It is well known how the mathematician’s interest in physical and astronomical 
problems have led to advances in mathematics itself. We cite here an example of a similar 
advance in mathematics which had as its origin a pure economic motivatiop. Not sur- 
prisingly, these new mathematical advances brought about have themselves stimulated 
the economic problem from which they sprang. 

The Menger seminar in Vienna on mathematical economics, amongst many other 
topics, concerned itself with that of the existence of equilibrium. This led to the paper 
by A. Wald [51]. Soon after von Neumann [42] in proving the existence of an equilibrium 
point for an economic system, found it necessary to extend the Brouwer fixed point 
theorem. Kakutani [28] then simplified von Neumann’s proof and cast the theorem in a 
somewhat different light. This lead to an even more powerful topological fixed point 
theorem by Eilenberg and Montgomery [16]. Begle [6] took up from there and gave a 
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very general fixed point theorem which subsumed that of Eilenberg and Montgomery. 
To complete the cycle, Arrow and Debreu [4, 12] have applied those theorems to obtain 
rigorous proofs of the existence of equilibrium points for fairly general economic systems. 

In conclusion, the author should like to offer an apology to those economists whose 
contributions to the field may have been overlooked in this paper. The limitations of 
space, as well as the author’s incomplete knowledge of the field, make such omissions 
unavoidable. An expository paper of this type cannot aim at complete coverage; rather, 
the author tried to convey an idea of the range and variety of mathematical results which 


have found applications in economics. 
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AN ITERATIVE SOLUTION TO THE EFFECTS OF CONCENTRATED LOADS 
APPLIED TO LONG RECTANGULAR BEAMS* 


BY 
C. A. M. GRAY 
The University of Sydney 


Summary. The analysis of the bending of thin deep rectangular beams by con- 
centrated loads has already been treated by L.N.G. Filon and Th. v. Kadrmdén. These 
authors consider an infinitely long beam and express the load as a Fourier integral, 
obtaining integral solutions for the stresses and displacements. In their integral form, 
these results are rather difficult to interpret, although F. Seewald, using Karmdan’s 
analysis, has calculated numerical values for the case of a single concentrated load. 

In this paper, a totally different approach to the problem is made. The loaded area 
is conformally transformed to a circle, and, using Muschelisvili’s transformation of the 
boundary conditions, the problem is solved in the plane of the circle. The solution is 
then obtained in the form of a complex power series. In this manner, a direct solution, 
to any required degree of accuracy, is readily obtained. 

1. Theory. Consider a long deep rectangular beam resting on two supports placed 
at the same level, and loaded by a weight placed midway between them. The weight of 
the beam will be neglected, and the beam will be assumed to be in a state of generalised 
plane stress, the plane of the mean stress being that of the length and depth. To determine 
the local effects of the concentrated loads on the stresses and deflexion it is convenient 
to move the supports to infinity and treat the beam as an infinitely long strip of depth 
2h as is shown in Fig. 1. This procedure is that adopted by L.N.G. Filon, H. Lamb, 
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z-plane C-plane 
Fra. 1. 


and Th. v. KArman [1, 2] in their solutions of the problem. 
Transform the strip ABCD to the unit circle, y, in the ¢ plane by the transformation 
— 9} < . 
po he 7, = —Aihe(t + 2/3 + 2/5 + ++). (1) 
7 aa 
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Following Muschelisvili [3], define two functions 2 (z) and w’(z) of the complex 
variable z by, 
oc. +o, = RY(2), o, — o, + 2ir,, = —${z2''(2) + w’’()}, (2) 
where the dashes denote differentiation with regard to z and the bars denote the con- 
jugate function. Then, the boundary conditions in the z plane can be expressed in the 
¢ plane by 
Q[f(o)] + flo) Qhe) + wf] = Fi + iF. , (3) 


where, 
z2=f() = Dour’, 


is the transformation relation, o is the value of ¢ on the unit circle, 





ag) = 


F, + iF, is the transform of f; + if. , f: + if. = 4ifo (X, + iY,) ds, X,, Y, are the 
components of the boundary tractions, applied in the z plane, and s is the distance 
measured along the boundary in the z plane. 

Since 2/(¢) is analytic within the unit circle, we can express it by the Taylor’s Series, 


o(¢) = >> br’. (4) 


Let us multiply both sides of (3) by (1/2 wi) do/(o — £) and integrate around y. 
Using Cauchy’s Theorem and, neglecting the constant, we have, 


. 


Q[f(¢)] = (1/2272) | (F, + iF.) do/(o — §) — (1/2z7) IF s(0) %@) da/(o — ¢) (5) 
Similarly from the relation conjugate to (3) we obtain, 
w’ [f(¢)] = (1/2z71) | (F, — iF.) da/(o — §) — (1/2n1) | f(a) Qa) da/(o — £).7Er (6) 


Substituting (4) into (5), carrying out the integration and noting that 
d | 


dé 2m. 


can be expressed as a power series, >-,.o p, ¢", we obtain the following infinite set of 


| (F, + iF.) = : 
te 


Y 


equations 


>. + v)u,.:50,-, + 1 +n) > Nisneste = De (7) 


Solving equation (7) gives 2! . To complete the solution, we require w’(z) which is 
established in exactly the same way, using (6). The formula for w’(z) is then 


. 


w'(2) = (1/2mi) | (F; — iF.) doo — ) — Der’, (8) 


where the c’s are given by the relations, 


Cr = >, Uda, (9) 


A full and complete account of the derivation of (7) and (8) is given in Reference 4. 
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2. Solution. Returning to the original problem, consider the given loading as a 
combination of the following two systems, where in system B, the points L,M,N,P are 
supposed to be very remote. 




















| v2 
D 
A = C 
x 
B 
[ve 
System A 
{74 |v2 {v4 
P D N 
L B M 
1/4 \i/2 4 
System B 
Fia. 2. 


Case A. Load function: 
along BCD, F, + iF, = —2i, 
along DAB, F, + iF, 


Il 
ad 


Then, 


(1/2nt) / (F, + iF 2) do/(o — {) = —x “log (on — $)/(on — 9) 
5 (1 /2ri) [ (F, + iF.) de/(o — t) ia 2 (1 + )-* 


(1 2771) [ (F, = iF.) da/(o am ) _ .- log we Fey iy/(er”? _ t). 
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Equations (7) become 


2b, + .33333b, + .2b, + .142857b, + .11111b, + --- = —(2h)"', 

2b, + 1.6b, + .42857b, + .33333b, + .27273b, + --- = (2h) 

2b, + 1.7108), + 1.5556b, + .45455b, + .38462b, + --- = —(2h)™', 

2b, + 1.778b, + 1.6364b, + 1.5385), + .46666b, + --- = (2h), ete. 


Using an iterative process to solve these equations we obtain for the 6th and 7th 
iterations the following values: 





Iteration bo Xh b, xh b, X h fe be Kh a 


6 — 32113 .91461 — 1.0132 .98457 — .98724 98556 
7 — 34898 88018 — 1.0381 .97263 — 1.01438 1.0143 











From these we obtain 


Q% = —.329/h + .895¢7/h — 1.025¢°/h + .979¢°/h + ---, 
a. a : (s 
=— 5 = + .171/h — .105¢°/h — .025¢°/h + ---, ) 
dQ’ d¢ lr - 2 2 ol | . 
27 = - = sct1l— £ = nz — 21 — .10¢ — .126¢ oe Pe 
dt dz th r(1 ‘ IF + on | : 263 a » 


Using these values in equations (8) and (9) we calculate w’(z) and w’’(z) as follows: 


l (= — ¢' tf 
(> = oO — = 
aia log 1—¢ 2 1 + 


? 21 a e ee ee 
3+ {— .55014¢ + .771667¢ 
$ 





T T 
+ .043¢° + ---], (c) 
, — ¢°| 1.2854 4 , — i = , 
w’'(z) = oA ! ¥, -£l a — .11763 + .1075¢? + .107¢* + | (d) 
Case B. Load function: 
For BM | MCN ND DP PAL LB 
F, + iF, | 0 | —1 —2i 0 —1 —2i 





Then, 


- [ (F, + iF _ eS £ log (ex = i)(s2 = f)(se — f)(ee a ‘)"] 
ant J, ie it ae e\ow —ft x —~— 5 Gp — S/\o, — ¢ 


Now when L and P move to A o, = op = 1, and when M and N move toC oy = oy = 
—1. Thus, 


y do —] —|] — ¢? 
ny . - _ a a ea ] » a = 
2Qri i (Fi iF) o—f¢ T 6 l-—-¢ 


Y 
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The constants b are purely imaginary and equations (7) become 


b, = 2b,/3 — 2b, /5 — 2b;/7 = 2b,/9 + oc = —I/h, 


b, + b; — 4b,/5 — 4b3/7 — 4b;/9 — 4b;/11 + --- = 0, 
b, + bs + bs — 6b,/7 — 6b;,/9 — 6b;,/11 — 6b,/13 + --- = —I/h, ete. 
Assume 
2? = irb,z/4h a f° a [Pw io —) -- 
i.e., assume that the b, are of the form 
b, = D,, 
bs = b,/3 + fs, 
b, = b,/5 + fs, ete. 
Then since 
2 2 2 2 
'-$~g5"Fs"ts° °° * > 
] 4 4 4 
de tue y tat i oad * tear 
1 l 6 6 6 
i++ sis" 38" in’ "> 


b, can be eliminated from the equations and we are left with the following system for 


the values of f: 


— .2f; — .142857f; — .111111f; — .090909f, + --- = —1/2h, 
107143f; — .111111f; — .090909f, — .076923f, + --- = 0, 
055556f; + .075758f; — .076923f; — .066667f, + --- = —1/6h, ete. 


To bring these into a form suitable for the iteration process subtract each equation 


from the one above it. Thus, 














— .307143f; — .031746f,; — .020202f, — .013986f, + --- = —1/2h, 
051587 f; — .186869f, — .013986f, — .010256f, + --- = 1/6h, 
.021365f, — .027681f, — .135256f, — .007843f, + --- = —1/6h, etc. 

From these we obtain 
fs + t/h 1.581 + .004 fo + t/h | —.72+ .01 
fs + t/h | —.581 + .001 fu + i/h | 1.25 + .01 
fr = ish 1.33 + .02 fis + t/h | —.77 + .02 
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Thus 
hQl = ibzw/4 — if¢/(1 + &°) + ig(1 + 5810? + .419¢* + .33¢° + ---) (e) 


and 


W202!’ = ibshe/4 + ; (i~ py/0 + 7 — 7 (1 — g*)[1 + 1.7432 


+ 2.095¢* + ---] (f) 
From formulae 8 and 9 and using the tabulated values of /, we have, 
w'(z) = ibe’r/Bh +07 log (1 t+ f)/1-&)4+7-04+ 2)" 
+ (4/m)[eo + .972¢° + .749¢° + -+-], (g) 
and 
hw’'(z) = i¢/(1 + &°) + (nr /2)¢(01 — )/0 + ©) + iben/4 
+ i¢(1 — ¢°)(1.944 + 2.996%? + ---). (h) 
The stress distribution for which Q!’ = ib,zr/4h, w’’(z) = ib,zr/4h is that due to 


pure flexure, provided b, = —7¢B,/h, where B, is real. From statics, we know that there is 


an infinite bending moment at the centre section due to the loads at infinity. But, from 
statics, at a section where the support loads are applied, the normal stress distribution 
must be equivalent to a moment of h/27 as shown in Fig. 3. As this moment is applied 


1/4 
IAT)IAT 


oo —— to 


lIAIZANIAT 


1/4 


Fia. 3. 





° ° —_— ° e ° . > 2 
effectively at infinity, it must be included in the b, term. For as ¢ — + 1, RQ} — B,z/4h’. 
Hence, B, is made up of two parts, that due to an infinite moment which will be repre- 
sented by b’ and that due to a moment of amount h/27. Thus 


B, = b! + 3/x’ 
and 
hQl = b’ew/4h + 32/4eh — i¢/(1 + ©) + ig + 581g? + 41964 +--+), (i) 
w'(2) = b’e'e/Bh° + 32*/8ah! + x * log (1+ )/1- o) +4 
—(1+ 9°)" + (4/)@ + 972° + ---) (i) 
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Equations a to f now give the solution to the problem. To calculate the stresses at 
any point we use equation (1) with the values of 2! ,.2/’ and w” obtained by substituting 
the corresponding values of ¢ in the relevant equation a to j. Fig. 4 gives the distribution 





Ty) 7’ 


*Y x20 































































































A 
Z 
, 4 xn[oy) 
-3 +2 . 2 s 4  xhlo,) ig | 
2 3 4 y/h 6 
VA : 
sie te a 
h 
nl hx ee 
LAL de 
} 
c= 
~ 
~ 
+-.8 
4-1, 





Fia. 4. 


of normal stresses across BD and the shear stress distribution along the y axis. 
The deflections u and v are readily obtained from formula (10) 


Su(u + iv) = (3 — )/(1 + )O® — 2Y@ — we, (10) 
where yu is the shear modulus; » is Poisson’s ratio, u, v are the components of displace- 


ment in the z and y direction respectively; and D is u + w. 

We are particularly concerned with the displacement and curvature of the y axis 
and for this case some simplification of 10 is possible. 

Differentiate both sides of (10) with regard to y, then, 


8u > = 10(3 — 9) (1 +») — 19'@ + icO"@) + B"®. 
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Since it is only the loading condition of case B that will have any effect on the centre 
line we have for x = 0, 





= = _ __ ¢2\2 1 ¥ __ 
iQ") + W'@) = sy: (i= £) +i. $47.08 


+ .944&/h — .024¢/h + 


and therefore 


Ou 4 —rb’y l E } 3 | us 

— = ens iS Se eres caw! tne A80¢ er en“ 
Se ou | 1+ | we TA Tee CORT OR +: ane” ¥ 

L= ay to ee eee 
(Lat +i Taeetont ‘ey. 5 + (/h)[.944 — .024¢ + ---]. 

To compare this result with the usual engineering form let u, be the deflection of 
x = 0 and, considering the case y > 0 let § = —t and FE denote Young’s Modulus. Then, 
du,_ rb’ y , 8 (uv) 1 l t 1 — 
ee se L =» (2) — — — , — 6961 — : 
de th RTS NA) BR BR Tee te OO — 


ai : ian 1 |e y (3 - ty t 
4800 — .405¢* + .358 --J-sol- ———, ———; 
+ .480¢ 405t° + .358¢° + ---] Suh E nh \iae) TT+? 


1-¢ 
Se ee ao 


———55 + .944¢ — 024° + --- 
; Ti + -944t — 0240? + | 


Using Euler’s summation process we have for the first series on the right-hand side 
2 
.696t — .608¢° .480¢° s+ = 696 s(t ) ° 
996 ) + .480¢° + 196 > + «08 xy + 


Then if A is the correction to the usual engineering deflection, 
dA —1 t (. (—! y 
— = ——.——_, 696 —— { ves 
dy ~ Eh 1+@7* e\? a rr eel Ce 


Suh k h 2 en 2 + 2 t (] J. ry + .944t .0242 oe 


Hence the correction to the usual engineering curvature is 





Pa DE as - ign + tel) + --) 
K= dy "2 E .696 1.216 5 rer + .224 ve a 
a me eS sy ri-? 1. 1-f _ wy 81-?) 
4vrE\l+¢ 32, bh? podilen (1+ ¢)’ 4h l 24 
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and, 
_-l, 1+, 1 y , 1.656 1 _ 088 _1 | 
4= 77 "ie :+4{ s70% es Tae 7 i¢+e@t 
_tly -». .236 ie ee 1.568 — 
Suh TeetE: =, og (1 — £) 8u rE 


Figure 5 gives a plot of K with y assuming vy = 0.3 and shows that for y > h, K is 
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negligible. Hence the engineering approximation, that the curvature is proportional to 
the bending moment, is quite exact at a distance from the load greater than half the 
depth of the beam. Making the assumption that for y > A the curvature of the centre 
line is proportional to the bending moment, we can then calculate the correction 4 which 


is also given in Fig. 5 
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MODIFIED STURM-LIOUVILLE SYSTEMS* 


BY 
W. F. BAUER 
University of Michigan 


1. Heat conduction problem. Consider the steady-state temperatures U(z,z) of a 
thin slab of material in the shape of a square with edges at z = 0, x = 1, z = 0, and 
z = 1. Let the edge at x = O be kept at temperature zero and each point of the edge 
at z = 0 kept at temperature F(x). Suppose the edge at z = 1 is in perfect thermal 
contact with a thin strip of anisotropic material the coefficient of conductivity of which 
is much greater in the x direction than in the z direction so that the temperatures in 
the strip can be regarded as a function of z alone. Suppose that an analogous situation 
holds for the edge at z = 1 so that the temperatures in that strip may be regarded as 
a function of x alone. If at the faces heat is transferred according to a linear law into 
an external medium at temperature zero and the coefficient of heat transfer r(x) is a 
continuous function, the boundary value problem may be written** 


U., + U.. — cr(x)U = 0, 
U(O, z) = 0, 

U,(1, 2) — k,U,.(1, 2) 

U(x, 0) = F(x), 

U(x, 1) — k,U,,(z, 1) = 0. 


In this statement of the problem 1/c is the conductivity of the slab, k, = cq, , kz = cq2 , 


where g,; and q, are the coefficients of conductivity of the strips along x = 1 andz = 1 


0, 


in the z and z directions respectively. 

The assumption made above of a narrow (relative to the size of the slab) strip of 
material whose temperature varies in only one direction may be closely realized in 
practice. For example, consider a material such as concrete which has imbedded in it 
pipes containing moving water. Also, consider a laminated material such as a brick wall. 

In the solution of the problem we are led to the classical method of separation of 
variables, since, because of the boundary conditions at x = 1 and z = 1 and the non- 
constant coefficient in the differential equation, there is little hope for success in the 
employment of a finite Fourier transform or the Laplace transform. In deriving solutions 
of the homogeneous conditions of the form U(z,z) = y(x)w(z) we obtain the following 
eigenvalue problem in y(z): 

y’"(x) + [A — er(x)]y(a) = 0, 


y(0) = 0, (2) 
~s ckyr(1)y(l) — y’(1) — ky’) = 0. 


*Received July 21, 1952. This paper is a condensation of parts of the author’s doctoral dissertation 
“Modified Sturm-Liouville problems and associated integral transforms’ written at the University of 
Michigan, 1951, under Professor R. V. Churchill whose advice and suggestions are gratefully acknowl- 





edged. 
**Throughout this paper function values at an end point of an interval shall mean the limit the fune- 


tion assumes as the point is approached from the interior of the interval. 
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In.order to justify the classical method of satisfying the non-homogeneous boundary 
condition U'(7,0) = F(x) of problem (1) by combining solution of the type y(x)w(z), 
an expansion theorem for the eigenvalue problem (2) is required. In this paper we shall 
deal with properties of and expansion theorems for problems of which type (2) is typical. 

2. The modified Sturm-Liouville problem. The eigenvalue problem (2) is a special 
case of the problem consisting of the Sturm-Liouville equation 


y’'(x) + [A + p(x) ly(x) = 0 (3) 


IIA 


in which p(x) is a real, continuous function for 0 S x 1, with the boundary conditions 


Uoly] = ay(O) + azy’(0) + asy’’(0) = 0, 
(4) 
U,[y] = by) + boy’) + by’) = 0 


in which the a’s and the b’s are real constants. In view of the differential equation (3), 


— 


the boundary conditions may be written in the alternative form 


Uyoly] = {a, — as[A + p(O)]}y(O) + a.y’(0) = 0, 
(5) 
Unily] = {b, — bsA + p()]}y(1) + bay’(1) = 0. 


In this form the boundary conditions are of the form of the usual Sturm-Liouville 
problem except for the appearance of the parameter ) in the coefficients. 
For two eigenfunctions y; and y; corresponding to distinct eigenvalues A, and }, 
we obtain by a familiar procedure the Green’s formula 
al 


(A; — Aj) yy; dx = [yys — ysyilo . 


0 


Supposing that a. ~ 0 and b, ¥ O we substitute from (5) into the right member and 
obtain the orthogonality relation 


al I 
| yy; de — os y,(0)y,(0) + - y.(1)y,(1) = 0. (6) 
‘ 9 Vo 

It is easily verified that if a, = 0 or b, = O the term of the orthogonality relation in- 


volving that quantity does not appear. 

Expansion theorems for various cases of the problem (3)(4), or the equivalent problem 
(3)(5), have been established by Langer [7], Gaskell [4], and Churchill [1]. Boundary 
value problems giving rise to the eigenvalue problem are described in the first two of 
these papers. Among these three papers, only Churchill’s deals with the problem* 
involving the differential equation with non-constant coefficients. In this paper we shall 
establish an expansion theorem for the problem having non-constant coefficients in 
the differential equation. The author believes that the theoretical development leading 
up to the expansion theorem as well as the theorem itself helps fill the gap in theory 
between eigenvalue problems which involve the parameter in the coefficients of the 
boundary conditions and those which do not. 

















*Another paper by Langer [8] presents, in another way, derivations of some properties of eigenvalue 


problems more general than ours. 
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3. Definitions and properties of normal and positive-definite eigenvalue problems. 
We make the following definitions: 

Definition 1. A function u(x) shall be called a “V-function”’ if it is real, not iden- 
tically zero, of class C’, and satisfies the boundary conditions (4). 

Definition 2. The inner product corresponding to the eigenvalue problem (3), 


(4) for two bounded integrable functions f(x) and g(x) is defined as follows: 
‘ a3 b; 
if, ol = [ fod — % fOg(0) + F f(A)g(1). 
Jo a2 2 


We note here that the inner product is defined for any two V-functions. Also, the ortho- 
gonality relation (6) is simply [y,; , y;] = 0 and the Fourier coefficient with respect to 
the function f(x) and corresponding to the eigenfunction y,(z) is written formally as 


[f, ys) 
[yi ’ yi] 


The inner product has the important properties of linearity and commutativity with 
respect to its functional arguments. 

Definition 3. The eigenvalue problem (3), (4) is “normal’’ if for every V-function 
u(x) it is true that 


[u, u] > O. 


It is easy to verify that a necessary and sufficient condition that the problem (3), (4) 
is normal is that a.a, S O and b,b, = 0. 

Theorem 1. If the eigenvalue problem (3), (4) is normal, it has only real and simple 
eigenvalues and real eigenfunctions. 

Proof. Suppose \ is a complex eigenvalue, y the corresponding eigenfunction, 
and the imaginary part of \ is not zero. Then ¥, the conjugate of y, is an eigenfunction 
and, according to the orthogonality relation (6), we have [y,7] = 0. But this is a contra- 
diction because the eigenvalue problem is normal (a,a; S 0, b.b, 2 0) and the integral 
term of the inner product is greater than zero while the last two terms are greater than 
or equal to zero. Therefore, if the problem (3), (4) is normal it has only real eigenvalues. 

Suppose next that u(a2) and v(x) are two eigenfunctions both corresponding to the 
same real eigenvalue. Since both of these functions are solutions of (3), the Wronskian 
of the two functions is a constant. If this constant is evaluated by means of either bound- 
ary condition in (5), it is found to be equal to zero. The functions u and v are, therefore, 
linearly dependent. That is, any eigenvalue of (3), (4) is simple in the sense that there 
cannot be two linearly independent eigenfunctions corresponding to it. It follows im- 
mediately that any complex eigenfunction u(x) + 7iv(x) is equal to (k + 7)v(x) for some 
k; any eigenfunction is real up to a constant factor. 

Let us now prove that the eigenvalues are simple in the sense that they are not 
double roots of the characteristic equation. Let y(z,A) be a solution of the differential 
equation (3) such that it satisfies the boundary conditions 


Uroly(x, )] = 0, y(1, A) = 1. 


The function y(z,A) and its first and second derivatives are continuous functions of x 
and analytic functions of \, for the system defining y(z,A) is an initial value problem 
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with a complex parameter \, and the coefficients of the problem are analytic functions 
of \. Therefore we may expand y(z,A) in the following Taylor series about the eigenvalue 
A; : 


y(z, +) = y(z) + D> A, (a, AA — AD’, (7) 
n=1 
where y/,;(x) = y(x,\,) is the eigenfunction corresponding to \; . Similar expressions can 


be written for y’(z,A) and y’’(z,A). 

In what follows we shall need to know that the coefficient of (A — \,) in the expansion 
of y’(x,A) is the same as A{(z,\,). This is true if 
0 0a 00 
——— > \) m —~— oe. ) 8 
an az YO ®) = az a YO ™ 8) 
which, in turn, is true if (0/0A)(0/dx)y(z,A) is a continuous function of x and X. Since 
y’(x,A) is analytic in X, 


a 4 - 1 F y'(x, r) " 
an YG) = oF I (@ — ve” 


where C is a finite closed path around \. This expression shows that because y’(x,A) is 
continuous in x and i, (z — A)’ is bounded from zero, and C is finite, (0/0A)y’(z,d) is 
continuous in x and ) and (8) is valid. In a similar way it can be shown that A{’(z,\,) 
is the same as the coefficient of (A — \,;) in the expansion of y’’(z,A). 
The problem which defines y(z,A) can be written as follows: 
y+ OA + py = —A— Ady, 
fa, — as[A; + p(0)]}y(O, ») + azy’(O, A) = as(A — A,)y(O, A), (9) 
y(1,») = 1. 
Into this problem we substitute the series (7) and the corresponding series for y’(z,A) 
and 7’’(z,A) and make use of the results of the preceding paragraph. The conditions in 
(9) hold for all \ and we may therefore equate coefficients of (A — \,) to obtain 
Ai’(x, Ai) + DAs + pa) JA, As) = —y.(2), 
a,A (0, Xj) — asy;(0), (10) 
A,(0, A,) = 0. 


We multiply the differential equation in (10) by y,(x) and the differential equation 
yi’ + (A; + p)y; = 0 by A,(z,A,), subtract the two, integrate and obtain 


1 


— | tye? de = [ty — Ail. 


Making use of the fact that the eigenfunction y; satisfies the boundary conditions 
Ux.0 = Ux. = 0 we suppose that a, ~ 0 and b, ¥ O and evaluate the right member 
by means of the boundary conditions in (5) and the two boundary conditions in (10). 
Thus we obtain 


tv vl — BP = 4 Oy aLAWG@, 20) (11) 
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Since the function y(z,) satisfies the first boundary condition in (5) for all \ the charac- 
teristic equation of the eigenvalue problem is 


Unly(x, )] = 0. 
The derivative with respect to \ of this left member when evaluated at \ = \, is 
U,,:[A.(2, \,)] — bsy,(1). 


We see from (11) that if this were zero we would have [y; , y;] = 0 which is a contra- 

diction in view of the fact that the eigenvalue problem is normal. The proof follows 

through with minor changes if a, = 0 or b, = 0. This completes the proof of Theorem 1. 
We next define the Rayleigh’s quotient R(u) for any V-function u(z): 


R(u) = — od + pu, u) 
lu, ul 


By taking the inner product of the left member of the equation y/’ + (A; + p)y; = 0 
with y; and solving for A, we see that R(y;) = A, ; the Rayleigh’s quotient for the eigen- 
function y,(x) (a V-function) is the corresponding eigenvalue \, . 

Definition 4. The problem (3), (4) shall be called “positive definite’’ if it is normal 
and if [u’’ + pu,u] < 0 for every V-function u(z). 

It follows that a positive-definite eigenvalue problem has only positive eigenvalues. 
By performing an integration by parts we find that if p(x) < 0 for 0 < x S 1 the con- 
ditions a,a, < 0, b.b, = 0, a,az S O, and b,b. = 0 guarantee that the problem is positive- 
definite. Other combinations of conditions which guarantee a positive-definite problem 
can be found. 

It is easily verified that the eigenvalue problem (2) arising from the boundary value 
problem (1) is positive-definite because in that problem r(x) > 0 for 0 S x S 1 since 
r is the thermal emissivity of the faces of the slab, and k, and k, are positive since q,, 
g2 and 1/c are positive since they are coefficients of conductivity of the materials in- 
volved. Hence everything that has been said and will be said concerning positive-definite 
problems holds for the problem (2). 

4. The Green’s function solution of the non-homogeneous problem. Consider the 
following non-homogeneous problem with the complex parameter X: 


y’'(x, +) + [A + p(x) ]y(2, »+) = f(a), 
U[y(a, 4)] = 0, (12) 
U,[y(x, A)] = 0, 


where f(x) isa bounded integrable function. Except for the function f(x) in the differential 
equation this is the eigenvalue problem (3)(4). We make the stipulation here that if 
in the boundary conditions a, = 0(b. = 0) then a; = 0(b; = 0). Making use of the 
differential equation we may rewrite the boundary conditions in the form 


Uyoly(x, )] = —asf(0), 
(13) 


Uy ly(x, A)] = —bsf()). 


It is well known (see, for example, [5] p. 257) that if \ is not an eigenvalue of (3)(4) 
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the solution of (12) with the boundary conditions in the form (13) is 


al 


y(z,») = | f@G(a, & d) d—é — asfO)Gd(z, ») — bs f(I)G,(z, 2). (14) 


In this equation G(z,£,\) is the Green’s function of the system and G,(z,A) is the solution 
of the corresponding homogeneous problem except that the right member of the first 
condition in (13) is replaced by 1. Similarly, G,(2,\) is the solution of the corresponding 
homogeneous problem except that the second condition in (13) is replaced by 1. 

Since the Green’s function satisfies the homogeneous problem we have 

Uyol[G(x, &, A)] = lim fa, — as[A + p(0)]G(a, &, dA) O=<£s I 
z—0 
—a,G,(x, —, )}= 0 

which gives 


s 


U,o[lim G(x, &, A)] = Uyollim G(a, & A)] — lim Ujo[G(a, &, r)] 
the last term being zero. Because G is continuous and G, has a unit jump across the 
line x = & we see from this last expression that 

U,o[lim G(x, —, \)] = a2. 

9 

Because of the properties of the Green’s function the function lim;.. G(7,£,A) satisfies 
the homogeneous differential equation corresponding to the one in (12) and the second 
boundary condition in (13). The function is also continuous and has a continuous de- 
rivative. Recalling the definition of G,(z,A) we see that, except for a multiplicative 
constant, this funetion must be identical with the function lim;.> G(v,é,A). More pre- 
cisely, if a, + 0 we have 

a . G(x, &, d) 

G,(z, A) = lim ; 


t— ? 
§-0 2 


In case b, ¥ 0 a similar argument yields 


G(x, g, A) 


G,(z, 4) = —lim 
: b. 


EI 
Therefore in view of (14)*and our definition of the inner product the solution of (12) 
may be written 
y(x, +) = [f(&), G(z, &, d)]. (15) 
According to Collatz [3], as long as \, is a simple eigenvalue the Green’s function 
has a simple pole at \ = A, and may be written 


Cry (x)y(&) 


G*(z, &, X) 
(X — d.) + G*(z, &, dr), 


G(z, —é, A) = 


where G* is analytic at \ = A, and C, is a constant. The quantity C,y,(r)y,(&) is, there- 
fore, the residue of G at the simple pole 4 = X; . It follows that 








lim {(A — AJ[y.6, G(x, —&, )]} = Cy (aly, y(é)]. (16) 


d 
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Using the fact that 
yi'(x) + [A + p(a)ly(x) = (A — Ady.(2) 


and the result (15) which gives an implicit solution for y;(x), we see that the left member 


in (16) is simply y,;(a). Therefore 


] 
C; = : 
[y(é), yi()] 
Thus the residue of G at \ = A; is completely determined and is found to be the term 
of the Fourier series of f(x) corresponding to the eigenvalue A; . If A,(¢ = 0, 1, --+ , n) 
is the set of eigenvalues inside the circle | \ | = R the meromorphic function G may be 


written 


eae - y (x)y.(€) 
G(z, =, \) = : : — Gk(a, &, Xd). 
2d (A — Ax)[yi(€), yi(E)] + Oe, 6 
where G# is an analytic function inside the circle | \ | = R. Also, in view of (15), the 
solution of the non-homogeneous problem (12) may be written 


; y(x)[f(, yl] _ ae 
y(a, A) = > A — rdly.@), Ol + [f(&), G(x, €, d)]. 


A note concerning an expansion theorem for V-functions. With the developments of 
this section and sections 2 and 3, the groundwork has been laid for the proofs of the 
existence of an infinite number of eigenvalues and an expansion theorem for V-functions 
for the problem (3)(4) which is positive definite. The proofs are analogous to those 
given by E. Kamke [6] for problems which do not involve the parameter in the boundary 
conditions. Because of the analogy with Kamke’s work and because the expansion 
theorem we shall prove in the next section is more general, the proofs are omitted. 

5. An expansion theorem by the Laplace transform. In this section we shall obtain 
an expansion theorem by the application of the Laplace transform to a problem in heat 
conduction. Let us consider the transient temperatures U(z,t) of a thin slab in the shape 
of a rectangular parallelopiped which has its upper and lower edges insulated. At the 
edge x = 1 the slab is in perfect thermal contact with a well stirred liquid the container 
of which is exposed to an external medium at temperature zero. The edge zx = 0 is 
kept at temperature zero. If the flat faces transfer heat into the surrounding medium 
according to a linear law and with thermal emissivity p(x), and if the initial temperature 
of the slab is F(x), the boundary value problem may be written as 


U(x, t) = U,A(2, ) > p(x) U(2, 0, 
U(O, 2) = 0, 
K,U(1, 2) + U.(, ) + K.U.(, }) = 0, 
U(x, 0) = F(x), 


where the units are adjusted so that the thermal diffusivity equals 1. In this statement 
of the problem K, = qA,/KA, and K, = hM/KA, and K, = hM/KA, where the 
constant g is the thermal emissivity at the end of the container, A is the coefficient of 
conductivity of the slab, h the specific heat of the liquid, M its mass, A, the area of the 
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slab in contact with the liquid, and A, the area of the container in contact with the 
external medium. We note therefore that K, > 0, K, > 0; also, p(x) > Ofor0 S x <1. 
Except for the variable coefficient of heat transfer of the slab which introduces the non- 
constant coefficient in the differential equation, this is essentially the problem con- 
sidered by Langer [7] and Gaskell [4]. 

The application of the Laplace transform 


| 


L{F(z)} = f(s = | e F(a) dx 


to the boundary value problem gives the following transformed problem: 
u’’(x, 8) — [s + p(x) ]u(z, s) 
u(O, s) = 0, (17) 


— F(x), 


(K, + Kys)u(1, s) + w’(1, s) = K.F(1). 


We note that the homogeneous problem corresponding to this is a case of the eigenvalue 
g 
problem (3)(4). From the results of Section 4 we see that the solution to this problem 
may be written 
u(x, s) = —[F(E), G(a, &, 8)], (18) 


where G(z,£,s) is the Green’s function for the system. From the definition of K, we 
noted that K, > O and hence the corresponding eigenvalue problem is normal; the 
problem has only real and simple eigenvalues s, . According to the results of Section 
4, u(z,s) has simple poles at those eigenvalues and the residues of those poles are the 
terms of the Fourier series corresponding to F(z). 

The existence of an infinite set of eigenvalues and an expansion theorem for the 
eigenvalue problem can now be proved with the aid of a close examination of the order 
properties of various terms involved in the solution (18). In this proof we shall deal 
with the Laplace inversion integral 
ga 

| e* f(z) dz 
1 


9 
8 a7 J y-iB 


L;*{f(s)} = lim 
and make use of the following theorem (see [2], p. 159): 

Theorem 5. Let v(s) be a function of a complex variable s and of the order O(s“) 
in a half plane ®(s) = yo where k > 1. Then the inversion integral L7"{v(s)} converges 
along a line y where y = yo to a continuous function V(t) which is independent of y 
and such that V(0) = 0. 

Let us now define functions u,(z,s) and u.(z,s) which are linearly independent solu- 
tions of the homogeneous differential equation corresponding to the one in (17) and 
are such that u,(0,s) = 0 and u{(0,s) = 1, while (K, + K,s)u.(1,s) + ui(1,s) = 0 and 


u,(1,s) = 1. We write the initial value problem which defines u,(z,s) in the form 


us’(x, 8) — su(z, s) = p(x)u.(z, 8), 


ux(1, 8) = —(K, + K;s). 
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By supposing that the right member of the differential equation is a known function, 
we find that the solution u.(z,s) may be written in the following implicit form: 





ee 8 =— los —s - K, + Ks 
u(x, 8) = sinh (4) + sinh [(1 x)s | corn (s’*) + “7 ] 
(19) 
l 2 1/2 
— an | sinh ((@ — Ds'*Ip@usl, 9) ae. 


It can now be shown that the function 
er 1/2 = on artes | 
M(s) = Max Soe sinh (7 ”_)un(z, #) exp [—(2 — 2s") | 
Oszsi | K, + Kys 
is bounded for R sufficiently large where s = Re‘’. Since M(s) is bounded we can see 
from (19) that 





2s'’* sinh (s'’*)u,(0, 8) exp (—2s'””) 
K, + Ks 


= [1 — exp (—2s'””)]{1 — exp (—2s'””) + [1 + exp (—2s8'”)]O(s''””)} + O(s"'””). (20) 


The zeroes of the function u,(0,s) are the eigenvalues as can be seen from the definition 
of u.(x,s). Since these eigenvalues are real we see from expression (20) that they cannot 
be large and positive. Therefore we set s = —p’ where p is real and find that the charac- 
teristic equation u,(0,s) = 0 takes the form 


sin ol sin p+ o(?) cos | + o(2) = 0. 


The applications of Rouche’s theorem incomplex variable theory gives the fact that 
there exists an infinite number of s,(n = 0, 1, --- ), that they are less than a fixed num- 
ber y) , and that when n gets large (—s,)'”” approaches nz. Since u,(0,s) is an analytic 
function it has only a finite number of zeroes in any finite region of the complex plane. 
Thus the existence of a denumerable set of eigenvalues is determined. 

Treating the function u,(z,s) in a manner similar to that above for u.(z,s), the fol- 
lowing result analogous to (20) can be obtained: 


Da 2 
K, + Ks 
In terms of the functions u,(z,s) and u,(z,s) the Green’s function of the system (17) 
appearing in expression (18) may be written 


((K, + K,s)u,(1, s) + wW(1, 3)] = 1 — ee"? + Of”) (21) 


. Fs sa —ulé, 8)ui(z, 8) < 
G(z, &, 8) u,(0, 8) ; egae<& 


om S)UAT, § 
ing seme i(E, 8)ua(z, 8) t<2z <1. 
u,(0, 8) 
If we next make use of the properties of the Green’s function and assume that F(z) is 
continuous and has a sectionally continuous first derivative, we may integrate expression 
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(18) by parts to obtain 
al 


su(z, s) — F(x) = F(1)G,(2, 1, s) — F(O)G,(2, 0, 8) — | G;(x, &, s)F(&) dé 


“0 


al r 
K.sF (lu 
= G(x, &, x) F(€)p(é) dé — =—+* 
" P(e) (Ko + Kis)u,(1, PS uN 8)" 
Because the quantities u.(0,s) and (K, + K,.s)u(1,s) + u{(1,s) appear in the denominators 
of the various factors and integrands above, the application of the results (20) and 
(21) yields 
. F(x) » 
u(x, 8) — = O(s ””*) Ris) 2=7,0<eS2 


Ss 


IIA 
| 


where ¢ is a fixed but arbitrary number. Thus we see that the conditions of Theorem 5 


are satisfied. Let us designate the residue of u(z,s) at s = s, by R,(x). If we suppose 
for the moment that the Laplace inversion integral of u(z,s) — F(a)/s can be repre- 


sented by the sum of the residues of the integrand we may write 


(22 


IV 


L;"{u(z, s) — F(x)/s} = z. e’''(R,(x)] — F(z), t 
and therefore, according to Theorem 5, 


Zz. R(x) = F(x) 


<0 
Because we have shown in Section 4 that the residues R,,(2) are the terms of the Fourier 
series corresponding to F(x), the proof of the expansion theorem rests on the proof 
that the equation (22) is valid, that the inversion integral can be represented by the 
sum of the residues of its integrand. 

Since the numbers (—s,)'”” approach nz as n gets large, the parabolas P,, given by 


| eee 
f= (n + 3) 2" csc 5 (gq = 1,2, +++) 


pass between the poles of u(z,s) when n is large. It can be shown that when s is a point 
on the parabola P, and is sueh that R(s) S vy 
F(x) 4 y 
u(x, s) — =| © . (23) 
8 (n + 1/2) 
for n sufficiently large where C is a positive constant; also, on the parabola P,, when 6 
is restricted to the range —7r <@S050< 7 


F(x) < D(6, , 
8 —~(n+ 1 


where the positive number D, as indicated, depends on 6, and ¢ and gets large as 0 
approaches x. In making these estimates we assume that F(x), in addition to being 
continuous and having a sectionally continuous first derivative, has a sectionally con- 
tinuous second integral around the closed path consisting of the parabola P,, and the 
line ®(s) = y equals the sum of the residues of the poles enclosed. The order properties 


u(x, 8) — 
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(23) and (24) are sufficient to guarantee that the integral of u(z,s) — F(x)/s over P, 
goes to zero uniformly in x for0 < ¢ S x S 1 — e with increasing n. In other words, 
the inversion integral of u(x,s) — F(x)/s may be represented by the sum of the residues 
of the integrand. 

Thus we have proved that F(x) may be expanded in the series of eigenfunctions 


of the problem 
y’'(x) — [A + p(x) ly(x) = 0, 
y(0) = 0, (25) 
(K, + K.)y(1) + y’(l) = 0, 


where A, and K, are constants with K, > 0, and p(x) is any continuous function on 


0 = x S 1. The expansion is 


F(x) = >> Avy(2), (26) 


i=0 
where 


_ LF), y(z)] 

~ [y(a), y(a)] 
in which the inner product is the one corresponding to problem (25). Our main results 
may be summed up with the following theorem: 

Theorem 6. The eigenvalues \, of problem (25) are real and simple, and there 
exists a real number yp such that A; S yo for all A; and the numbers (—),)'” approach 
ix as i gets large. For any continuous function F(x) such that F’(x) and F’’(x) are sec- 
tionally continuous on 0 S x S 1 the expansion (26) is valid on 0 < x < 1 and the 
convergence is uniform on any closed subinterval of 0 < x < 1. 

If in (25) K, = 0 the problem reduces to a standard type Sturm-Liouville problem 
for which expansion theorems are well known (see, for example, [2], p. 259). Although 
(25) is a special case of problem (3)(4), only minor changes are necessary in the proof 


given in this section to extend it to the more general problem. 
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APPROXIMATE SOLUTION OF AN INITIAL VALUE PROBLEM BY 
GENERALIZED CARDINAL SERIES* 


BY 
H. D. BRUNK 
Shell Oil Company, The Rice Institute, and Sandia Corporation 


1. Introduction. A problem which arises in many different contexts is to approximate 
one of a certain class of solutions, u(z,y), of a partial differential equation or integral 
equation by a function which interpolates its values on a line. More generally, let V be 
a vector space, whose elements are functions of a point P in a space X, so that the sum 
of two functions in V belongs to V, as does the product of a function in V by a constant. 
Consider the problem of obtaining a function of V assuming given values a, at points 
P, , for k belonging to a set, J, of integers. If a family {A,(P)} of functions of V can be 
determined so that A,(P,) = 1, A,(P;) = 0 for 7 # k (j,k e I), then such an interpolatory 
function is given formally by the sum } a,A,(P). Such functions A,(P) can some- 
times be defined as follows. Let, for each point P, g(P) denote an element in a Hilbert 
space H. Let (a,8) denote the inner product of two elements, a and 8, of H, and let 
there exist elements v, of H such that (¢(P,),%) = 1, (o(P;),v;) = 0,7 4 k, 7k e I. 
Then one may put A,(P) = (¢(P),v,), if this inner product belongs to V. Suppose, for 
example, that as functions of t, g(t;P) and »,(t) belong to the class L, (have integrable 
square) on an interval (a,b) (-2 <a <b < o), that the integral f° o(t;P;)v,(t) dé 
is 1 forj = k and Oforj ¥ k, j,k e I, and that the integral belongs to V. Then one may 
put A,(P) = f° o(t;P)v,(é)dt. If, as functions of t, the functions {¢(t;P,)} form an ortho- 
normal system over (a,b), one may set v,(t) = o(t;P,). 

We shall consider only the case where P represents a point in the zy-plane of a set E 
containing the x-axis, where »v,(t) = exp (ckt), and where $(¢;z,y) is defined for (z,y) e E 


so that 


g(t; x, 0) = exp (—iz?). (1.1) 
We define 
Arle.) = 3 . . o(t/\s 2,9) exp Gk dt” (1.2) 
for each integer k, each positive number X, and each point (z,y) in E. The series 
YE [OAs a) (1.3) 
becomes, for y = 0, 
PA (kX) sin k (x — a | f : (x — kA), (1.4) 
which is the cardinal series associated with values f(kA) at points x = kA (k = ---, 


—2, —1,0, 1, 2, ---). If a function U(z,y) is given for y = 0, U(xz,0) = f(z), the series 
(1.3) gives, formally, a function defined on E coinciding with f(x) at points x = ka 
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(k = --- , —2, —1, 0, 1, 2, ---) on the z-axis. The advantages of the cardinal series 
as an interpolatory function are well-known. The cardinal series (1.4), under suitable 
conditions, yields a ‘‘smooth” function of z, in the sense that its Fourier transform 
vanishes outside the interval (—2/d, 2/A) ({1]; ef. also [2]). There is also a ‘‘consistency”’ 
property: a new cardinal series associated with values of the cardinal series (1.4) at 
equally-spaced points having adjacent points not farther than \ units apart coincides 
with the cardinal series (1.4) ([3], [4], also [2]; for further references to the literature 
on cardinal series see [2], also [5]). In 1908 de la Vallée Poussin proved approximation 
theorems [6] giving conditions sufficient in order that the series (1.4) will approach 
f(x) as \ — 0. It is our purpose here also to obtain an approximation theorem: to obtain 
conditions sufficient in order that the series (1.3) should converge, and should approach 
U(x,y) as \ — 0, where U’(z,y) is a function of V, defined on F, uniquely determined by 
its values f(x) on the z-axis. In §3, the method and the theorem of §2 are applied to 
the problem of approximating the temperature at a certain instant in an infinite insu- 
lated rod in terms of its temperatures at a later instant. As is well known, this problem 
has also a probabilistic interpretation. If a random variable z having a known frequency 
function is known to be the sum of independent random variables, x and y, the random 
variable y having a normal frequency function with mean 0 and standard deviation co. 
§3 applies, and yields a method of numerical approximation to the frequency function 
of x. In order for the method to apply formally, it is not necessary that y be normally 
distributed. 

2. Anapproximation theorem. In this section standard techniques of Fourier analysis 
are used to prove an approximation theorem, Theorem 2.1. 

We observe that the function sin [r(x — kd)/A]/[x(x — kA)/A], and the function 
which vanishes outside [—7/\, 2/A] and which is given by \ exp (ikAt) on [—7/A, 2/Al, 


are a Fourier transform pair; i.e., 


— are = os |. : exp (tkAt) exp (—zzxt) dt. (2.1) 
If the series 
r 2 f(kX) exp (tkX2) (2.2) 


converges uniformly, then it is clear from (2.1) that the series (1.4) converges (through- 
out this paper, integrals and series extending from — © to © are said to converge if 
the Cauchy principal value exists, and converge uniformly if >>”, or f”, converges 
uniformly as p — ©). It is desirable, therefore, to determine conditions on f(x) implying 
that the series (2.2) converges uniformly. It is simpler to state such conditions as con- 
ditions on the transform, F(t), of f(x). We have the following lemma. 


Lemma 2.1: If 


F(t) is continuous and of bounded variation on (—@, @), (2.3) 
if, for each positive number i, 


>-F._. F(t + 2mj/d) converges uniformly and absolutely for | t | < /d to a func- 
tion F,(t), (2.4) 
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then the integral 1/2 J“. F(t) exp (—iazxt)dt converges to a function f(x), and the series 


2 > f(D) exp (ikAd) (2.2) 


converges uniformly to F(t); moreover, for | t | < m/d, 
lim F,(t) = F(d, (2.5) 
as 1/X — © through integral values. 


Proof: We have 


a(2p+1)"/A P w/r 
N | F(t) exp (—ikAt) dt = x» >> [ : F(t + 2mj/X) exp (—tkat) dt 


2p+1)a/d j=-p ¥—ar/ 


[ 3 Flu/d + 2xj/d) exp (—iku) du. 


vJ—r j=-p 


Il 


By hypothesis, >>?__.. F(u/A + 2xj/d) converges uniformly and absolutely to F,(u/A) 
for |u| < x. Since d is arbitrary (positive), it follows that (1/27) f°. F(t) exp (—izt)dt 
converges to a function f(x), and that Af(kA) = (1/2r) f", Fy(u/d) exp (—iku)du. 
We observe that the total variation of F,(u/A) on the interval | u | < 7 is not greater 
than the total variation of F(t) on (— ©, ©). Moreover, F,(u/A) is continuous (and of 
period 27, as a function of u). Hence its Fourier series, > oa S(kA) exp (iku), con- 
verges uniformly ({7], p. 42) to F\(u/A), so that \ >of__. (kA) exp (ikAt) converges 
uniformly to F(t). Hypothesis (2.4) clearly implies that, for | | < 2/A, lim 5 F,(t) = 
F(t); for let ¢ be fixed. We have | F,(t) — F(t) | < doy) | F(t + 2mj/d) | < DYosisy 

F(t + 227) | if \ < 1/N. But this latter sum is arbitrarily small for sufficiently large 
N. The proof of the lemma is complete. 

The hypotheses on F(t) of continuity and finite total variation may be replaced 
by others which imply the uniform convergence on | u | < z of the Fourier series of 
F,(u/d). This lemma is essentially equivalent to Poisson’s summation formula ([8], 
p. 33, ff.). 

We recall that 

g(t; x, 0) = exp (—iz4), (1.1) 


and 
ee ; 
Aralt, y) = 5 | o(t/A; 2, y) exp (ike) dt. (1.2) 
Using lemma 2.1, we obtain the following theorem. 
THEOREM 2.1: Let F(t) satisfy the following conditions: 
F(t) is continuous and of bounded variation on (— ©, ©); (2.3) 
>>*._.. F(t + 2nj/d) converges uniformly and absolutely for 
| ¢| < w/d to a funetion F(t); (2.4) 
to each point (x,y) in E corresponds a function of t, L(t;x,y), integrable with respect to t 
over (— ©, ~), such that, for each X > 0, 


Fy ()e(t;2,y) | < L(t;x,y), for |t| < w/d. (2.6) 


| 
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Then the integral (1/2r) f=. F(t) exp (—izxt)dt converges to a function f(x). The series 
>. S(RA)Az.a(z,y) converges to the function f,(a,y) = (1/2m) S72, Fi(t)e(tz,y)dt for 
(x,y) « E; fy(x,0) is the cardinal series associated with values f(kd), so that in particular, 
fx(hA,0) = f(RA) (F = +--+ , —2, —1,0, 1, 2, ---). Moreover, as 1/X © through integral 
values, f,(x,y) — f(x,y) for (x,y) in E, where f(x,y) = (1/2r) f°. F(He(t;z,y) at, a function 
coinciding with f(x) on the x-axis.* 

Proof: That the integral (1/27) f°. F(t) exp (—izxt)dt converges, is guaranteed by 
lemma 2.1. Also, by lemma 2.1, the series \ }>¢--. f(kA) exp (ikAt) converges uniformly 
for | t| < 2/d to F,(é). We have 


n ] aT n : 
> f(kA)A:z.(z, y) = on >> f(kd) exp (tkt)e(t/d; x, y) at, 
k=—n v¥—a7 k=—n 
X a«e/r n : 
=> | > f(kA) exp (tkduo(u; x, y) du. 
2r #—r/K k=—n 
Hence 
> f(R) Azz, y) converges to f(x, y), (2.7) 
where 
1 a@w/k 
fiz, y) = _ | _ Fy(be(t; x, y) dt. (2.8) 
By (1.4), f,(z,0) is the cardinal series assuming values f(kA) at pointsz = kA (k = ---, 
—2, —1, 0, 1, 2, ---). By (2.5), hypothesis (2.6), and Lebesgue’s bounded convergence 


theorem, we have 


ae sane x 
lim fi(z, y) = on # F(t)e(t; x, y) dt = f(a, y), 


which completes the proof of the theorem. 
Essentially, Theorem 2.1 gives conditions sufficient in order that an integral of the 


form {*.. F(é)g(i)dt = f2. o(t)dt f*. f(x) exp (ixt)dx can be approximated by replacing 
f(x) by the cardinal series associated with its values at points = kA (k = --- , —2, 
—1, 0, 1, 2, ---), \ > 0. The parameters (z,y) are mentioned explicitly in the above 


formulation with a view to the application in which it is desired to continue into a set 
E of the zy-plane a function given on the z-axis, which belongs to a certain class. 

3. Temperature on an infinite insulated rod. Let a be a positive constant, and let 
U(z,y) denote the temperature at the point with coordinate z on an infinite, insulated 
rod, at time y > —a. It is known that if U(x, —a) is piecewise continuous, bounded, 


*Note added in proof: The hypothesis in Theorems 2.1 and 3.1 that F(¢) is continuous and of bounded 
variation on (— ©, ~) serves only to justify the term-by-term integration of the product of ¢(t/A; z, y) 
by the Fourier series of F,(¢). Accordingly, hypothesis (2.3) in Theorem 2.1 may be replaced by the weaker 
hypothesis that F,(¢) is integrable over [—2/d, x/X], if the further hypothesis that ¢(t/A; z, y] is of bounded 
variation as a function of ¢ (z, y, \ being fixed) on [—7, x] is added. Since these hypotheses are satisfied 
in the special case to which Theorem 3.1 applies, in that theorem the hypotheses of continuity and 
bounded variation on F(t) may be omitted. Correspondingly, in Corollary 3.1 the hypothesis that x“f(z) 
is absolutely integrable over (— ©, ~) for some a > 1 may be replaced by the hypothesis that f(z) is 


absolutely integrable over (— ©, ~). 
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and satisfies a Lipschitz condition, then the function U(z,y) given by 


i) 


Ue) = aap | Ue —a) em [-€ - N/A + a) de 


is the unique solution of the heat equation for y > —a which, together with its partial 
derivative with respect to x is bounded (| U(z,y) | <M, —-7 <z4< ~,y > —a, 
| Uayy) | < M,,|x| > 2,,y > —a) and which approaches U(x, , —a) at a point of 
continuity as (x,y) approaches (2) , —a) from above (y > —a). 

Let V be the class of functions u(z,y), defined for y > —a, such that 


ult, = ay pops | ue —O exp {-E- 9) AY +a} a G1) 
for y > —a. We take for the set E the set of all points (x,y) with y > —a. Let us suppose 
that U(x,y) e V, that U(z,0) = f(x) is known, that U(x, —a) is absolutely integrable 
over (—, ~), and that it is desired to express U(z,y) for y > —a in terms of f(z). 
We observe that for fixed y > —a, U(z,y) is the convolution of the functions U(x, —a) 
and exp {—2"/4(y + a)}/2[x(y + a)]'”. Thus if V(t,y) = f2. U(z,y) exp (ixt)dx, we 


have 


V(t, y) = Vit, —a) exp {-—(y + a)} for y>—a. (3.2) 
Hence 
V(t, 0) = V(t, —a) exp (—at’), (3.3) 
and 
V(t, y) = Vit, 0) exp (—y?’) for y> —a. (3.4) 
Since U(z,y) = (1/2) [2. V(t,y) exp (—itx)dt, we have 
U(x, ») = & [exp (—ite — yt) at [ UG, 0) exp (igh) de (3.5) 
iw Je Jaw 
for y > —a. This formula gives U(z,y) in terms of U(z,0) = f(x), but in a form un- 


suitable for numerical approximation. If y > 0, we may interchange the order of inte- 
gration in (3.5), or use (3.4) directly, to obtain 


Ula, ) = acy | UG, 0 exp {—@ — 2)*/4y} a, (3.6) 


but this formula is not available for y < 0. The method developed in §1 and §2 applies, 
however, and yields an interpolation formula, 


fia) = D> $0) Arale, W, (3.7) 


which, under suitable conditions, approximates U(z,y) for y > —a. 
Let 


g(t; x, y) = exp (—itz — yt), (3.8) 


and 


Anale, ) = 5- | olt/d; 2, 9) exp (it) at (1.2) 
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We have 
A,,,(2, y) = ®(2/vA — k, y/d’), (3.9) 


where 


l . 
P(x, y) = on | exp (—ilz — yt) dt. (3.10) 


Applying Theorem 2.1, we obtain the following theorem. 


THEOREM 3.1: Let U (x,y) be a function of the class V; i.e., let 


< ] — . , 
U(z,y) =s— 3 | U(é, —a) exp {-—(€ — x)*/4(y + a)} dé for y> —a. 


2[r(y + a)]' 


Let U(x,—a) be absolutely integrable over (— ©, ~). Let F(t) = V(t,0) = f°. f(x) exp 
(txt)dt, where f(x) = U(x,0), and suppose that F(t) is continuous, and of bounded variation 
on (—o, ~). Then the series > ae f(kA)®(x/X —k,y/d*) converges to a function f(x,y) 
of V such that f,(x,0) is the cardinal series associated with values f(kd), so that in partic- 
ular f,(kX,0) = f(kd) (k = «++ , —2, —1, 0, 1, 2, ---). Moreover, as 1/X — © through 
integral values, f,(x,y) approaches U(x,y) (y > —a). 


Thus, if a sufficiently small unit interval is chosen on the x-axis, the formula U (x,y) ~ 
fi(z,y) = } ea . f(k)®(x — k,y) provides a method of numerical integration of (3.5), 
for negative as well as for positive values of y. 

Proof: The above discussion, leading to equations (3.2) and (3.3), shows that 
F(t) = V(t,0) exists. Since U(z,—a) is absolutely integrable over (— ©, @), the func- 
tion V(t,—a) is bounded. Equation (3.3) then insures that the series . ae F(t + 2m7/X) 
converges unifirmly and absolutely for | ¢ | < 7/A to a function F,(t) which is O(exp 
(—at’)). Hence, for y > —a, Fy(t)¢(t;z,y) = O(exp {—(y + a)t’}). Thus the hypotheses 
of Theorem 2.1 are satisfied, the set EH being the set of all points (z,y) for which 
y > —a. We conclude that as \ — 0, f,(z,y) — f(x,y) = (1/2r) fF. F(te(tz,y)dt; 
but by (3.5) this is just U(2,y). It remains to show that f,(z,y) e V. One verifies easily 
that 


l ™ . 2 ‘ 
g(t; z, y) = Qnty + a}? | AGE —a) exp {[—(€ — a) /4(y + a@)} dé (8.11) 


for y > —a, i.e., d(t;z,y) e V for each ¢t. By Theorem 2.1, and (3.11), 


fi(z, y) = on | - Fy(t)e(t; x, y) dt 
= = ; 2 
~ 4r(y + a)” / _, oP (at) F(t) dt | , oP {—7t — (€ — x)" /4(y + a)} dé 
r ; SS is 
= aay + a)” in exp {—( — 2)" /4(y+q@)} axl A i Fy (e(t; &. —a) at| 


for y > —a, the interchange of integrals being justified by virtue of the uniformity 
with respect to ¢ of the convergence of the integral 


i) 


| exp {—7tt — (& — x)’/4(y + a)} dé. 
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But 
1 ew/h . 
fg —a) = 5- | Felts8, —a) at, 
Po 
hence 
iu. - l = @ aa a a ae ee 
fe) = aap | AG —a exp {-E- a)/My +a) dE (Y> -a), 


i.e., f,(2,y) e V. This completes the proof of the theorem. 
It seems desirable to state conditions on f alone, so far as is possible, sufficient for 
the conclusion of Theorem 3.1. To this end we prove the following lemma. 


Lemma 3.1: Jf x*f(x) is absolutely integrable over (—~, ~), a > 0, and if f*. 
U((2a)'*v + 2,—a) exp (—v*/2)dv converges uniformly with respect to x (in particular 
if U(x,—a) is bounded), then x*U(x,—a) is absolutely integrable over (—@, ©). 


Proof: Suppose the contrary; then if U,. and U_ denote the positive and negative 
parts of U(z,—a) respectively, either x*U,(x,—a) fails to be integrable over (0,~), 
or over (— ©,0), or x*U_(x,—a) so fails. Suppose the first eventuality occurs. We have 


ia" s U.(é, —a) exp {—(€ — x)*/4a} dé 


{.(2) = 


re) 


] ry 1/2 > ame )? ) 
ame | _ U.((2a}'"» + 2, —a) exp (—v"/2) do, 


aM ~ 


1 M+(2a)'/*e 
| z\* f.(2)dz = (2n)73 [ exp (—v’/2) dv [ | t — (2a)'”*v |* U.(t, —a)dt 


“N+(2a)*/?0 


1 S78 M 


1/(2a) 
a i ti al _ erry - 
> (Qn)'” i exp ( v, 2) dv [ ig 1) 7 A¢. a) dt 


If N > 0. Now (¢ — 1)* = #(1 — 1/0)* > (1/2*)t* on [N + 1, M] if N > 1, hence 


aM M 
| a*f.(z) dx > = Erf {1/(2a)'”*] [ t°U,(t, —a) dt, where 
Jy Jwet 
Erf(z) = 1/(2r)'” fi exp (—v?/2)dv. Since f¥,, t°U.(t,—a)dt is arbitrarily large for 
proper choice of N, M, x*f,(zx) is not integrable, contrary to the hypothesis. We obtain 
similarly a contradiction if z*U, is not integrable over (— ©,0), or if z“U_ is not inte- 
grable over (0,~), or over (— ©,0). 
In connection with the hypotheses of Theorem 3.1, we now make the following 


observations: 


(A) If f°. Ul(2a)'?v + 2,—a] exp (—v’/2)dv converges uniformly with respect to x 
(in particular, if U(x,—a) is bounded), and if f(x) = U(2x,0) ts absolutely integrable, so 
also is U(x,—a). This is an immediate consequence of Lemma 3.1, with a = 0. 

(B) If f(x) is absolutely integrable over (— ©, @), then F(t) = f*. f(x) exp (izt)dx 
is continuous. The proof is immediate. 

(C) If U(z,—a) >O2._. exp {—(anr/z)*} is absolutely integrable, then 7(t) is of 
bounded variation on (— ©, ©). 
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Proof: By (3.3), we have F(t) = V(t,—a) exp (—at’), where V(t,—a) = f2.. U(x,—a) 
exp (irt)dx. Hence 


N N 
> | Ft’) — FAt) | = dX | Vth’, —a@) exp (—at}’””) — V(t}, —a) exp (—at}’) |, 
i=M i=M 

where the subscript c (for cosine transform) denotes the real part of the corresponding 


function. This latter sum is given by 


i U(x, —a)[exp (—at}’*) cos xt)’ — exp (—at)*) cos xt}] dx J 
| 


i=M | J-@ 


which is not greater than 


a N 
| U(x, —a) | >> | exp (—at}’”) cos xti’ — exp (—at}’) cos xt} | dz. 
J-o i=M 


Now the total variation of exp (—aé”) cos at on (— , ©) is given by 2 }ox__.. exp 
{—a(nx/x)"}, hence the right hand member is not greater in absolute value than K, 
where K = {"°.. | U(x,—a) } aad » exp {—a(nr/x)*}dx. Thus F(t) is of bounded vari- 
ation on (— ©, ©). For F,(¢) we replace cos xt by sin zt. The total variation of exp (—at’) 
sin zt on (— ©, ~) is given by 2 )>~__.. exp { —al[(n + 4)x/z]"}. But bbe exp {—al(n + 
1)x/z]}"} < 02. exp {—a(nm/zx)*}, and >072_.. exp {—a[(n + 4)x/z}?’} < DUt--. exp 
{—al(n + 4)x/z]’} = 1+ Do@, exp {—al(n + 4)x/z]P} < 1+ > 2-0 exp {—a(nx/zx)’}. 
Hence U(x,—a) >ox--. exp {—al(n + 4)x/z]"} is also absolutely integrable over 


(— ©, ~), so that F,(¢) is of bounded variation on (— ”, @) also. 


(D) If there exists a number a > 1 such that x* f(x) is absolutely integrable over (— ~, ~), 
and if {°.. U((2a)'v + x,—a) exp (—v"/2)dv converges uniformly with respect to x (in 
particular, if U(x,—a) is bounded), then F(t) is of bounded variation on (— ©, ©). 


Proof: By lemma 3.1, x“U(z,—a) is absolutely integrable over (—”, ~). But 
>2._. exp {—a(nr/x)*} = o(x*) asx — @, for a > 1. To see this, we have only to 
observe that max, x’ exp {—a(nx/zx)*} = exp (—b/2) (b/2am*)’”(1/n)’. Hence x~° 
2. exp {—a(nr/x)”} < «> + 2 exp (—b/2) (b/2an’)’” D0%., 1/n”. On choosing b 


between 1 and a, it becomes clear that 27° )o°.-. exp {—a(nx/zx)*} ~O0 asx — o. 
Thus the hypotheses of (C) are satisfied. 
These observations yield the following corrollary of Theorem 3.1: 


Corouiary 3.1: If U(z,y) ¢ V, 2.e., if 


U(a, y) = ang a) |” I. Ug, —a) exp {-(E—2)*/4y+a)} for y>-—a, 
if U(x,—a) is bounded, and if x*f(x) is absolutely integrable over (—©, ©) for some 
number a > 1, where f(x) = U(z,0), then the series >o7--. f(kA)®(x/X —k,y/d*) con- 
verges to a function f,(x,y) of V such that f,(x,0) is the cardinal series associated with 
values f(kX) (k = «++ , —2, —1,0, 1, 2, ---). Mbreover, f,(x,y) — U(2,y) (y > —a), as 
1/A — © through integral values. 
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4. A probabilistic interpretation. It is convenient to perform a translation of the 
xy-plane parallel to the y-axis, so that the point (0,—y.) (0 < yo < a) becomes the new 
origin. Let f*(x) denote the function U(x,—yo), and h*(x) the function U(z,0). Formula 
(3.4) yields 


‘ wa : 
U(z, 0) = Wry.) fg U(E, —Yyo) exp {—(E — 2)" /4yo} dé, or 
Re a 
h*(z) = >waz | S*© exp [— — 2)*/Ayo} db. (4.1) 
47) 6 J_@ 


If h* is regarded as known, and f* as unknown, §3 gives a method of approximate solution 
of the integral equation (4.1). 
Put ¢ = (2y)'”’; we have 


h*(z) = - _—_ | f*) exp {—(é = 2)" /20°} dé, 

(29) o J ow 

so that h*(x) is the expected value of a function f*(w) expressed in terms of the mean, 
x, of the normally distributed random variable w with standard deviation o. Thus §3 
provides a means of approximating a function f(z), if the expected value of the function 
f(w) is known in terms of the mean of the normally distributed random variable w 
having known standard deviation o. 

If f(x) is a probability frequency function, then h*(x) is the probability frequency 
function of the sum of two independent random variables, one of which has the frequency 
function f(x), while the other is normally distributed with mean 0 and standard devia- 
tion o. Thus if it is known that a random variable z having a known distribution is 
the sum of two independent random variables, x, and y, y being normally distributed 
with mean 0 and standard deviation o, the method of §3 can be used to approximate 
the frequency function of x. 

It is not necessary that y be normally distributed in order to apply the method 
formally. Let x have the unknown frequency function f*(x), y the known frequency 
function g*(y); let x and y be independent, and let z = x + y have the known frequency 
function h*(z). Then h*(z) = f2. f*(x)g*(z — x)dx. Let H(t) denote the transform of 
h*, F(t) the transform of f*, and G(¢t) the transform of g*. Formally, we have H(t) = 
F()G(t), F(t) = H(t)/G(b), and f*(x) = (1/27) f2. exp (—izt)/G(dt f*. h*(€) exp 
(itt)dé. Now replace h*(£) by the cardinal series associated with its values at points 
k\(\ > 0, & integer): 


h(¢) = me h(ky) sin 5 (2 — a) /* (a — ka). 


k=—@ 


We obtain 


3 . */()- l \ af 1 - . / - ;]- 

falz) = DWBA) 9 a Guay °XP (—ixt/) exp (kt) dt. (4.1) 
Thus if we put ¢(t;z,y) = exp (—izt) for y = 0, o(t;z,y) = exp (—izt)/G(® for y = yo, 
and let E be the set consisting of the lines y = 0, y = yo , then equation (4.1) is the 
result of a formal application of the method to which Theorem 2.1 applies. However, 
it may be difficult, in individual cases, to verify the hypotheses of the theorem. 
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AN ASYMMETRICAL FINITE DIFFERENCE NETWORK* 


BY 
R. H. MACNEAL 


California Institute of Technology 


Introduction. Finite difference techniques have been used extensively in recent years 
in the solution of two-dimensional second order boundary value problems that have 
proved to be intractable by other methods. The differential equation is replaced by a 
system of linear algebraic equations, the solution of which gives the values of the wanted 
function at a finite number of points lying at the intersections of a gridwork. The use of 
regular polygons, either squares or equilateral triangles, in the formation of these grid- 
works has the desirable property that the equations associated with each node (inter- 
section) point have a particularly simple, symmetrical form that is identical for all 
interior points. There are, however, two troublesome problems which arise in connection 
with the use of regular polygons. The first of these arises when the region has curved 
boundaries. In such cases some node points near the boundary will be connected to the 
boundary by gridwork elements of irregular lengths, necessitating the use of special 
equations for these points. The second problem concerns the change of mesh size at 
points within the boundary. It is frequently uneconomical from the point of view of the 
labor of computation to use the same mesh size at all points. In the neighborhood of a 
sharp corner or near other types of singularities, the mesh size must be reduced if an 
accuracy is to be obtained that is comparable with the accuracy of the solution in parts 
of the region where the behaviour of the wanted function is more uniform. Both of these 
problems have received attention from writers on relaxation methods and it is with these 
problems that the present paper is principally concerned. A method will be described 
by means of which the coefficients of the system of algebraic equations can be computed 
for an arbitrary distribution of node points. The positions of these node points can then 
be chosen to fit the boundary conditions and other special requirements of each problem. 

In the construction of a finite difference gridwork to be used in the solution of physical 
problems, it is helpful to associate physical properties with the elements of the gridwork. 
Southwell and his co-workers have regarded the gridwork as a network of tensioned 
strings [1] while others have regarded the gridwork as a network of electrical elements 
{2}. In the case of a second order boundary value problem, this clear physical picture of 
the gridwork is lost if the differential equation is replaced by difference equations in- 
volving difference operators higher than the second order. In order to preserve the 
physical picture and to simplify the calculations, the higher order difference terms are 
usually regarded as corrections which are added in the final stage of calculation, if the 
relaxation method is to be used [3]. 

Another important reason for eliminating higher order difference operators arises 
in connection with analog computing devices, in which the physical picture of the net- 
work is realized. Electrical circuits have been used extensively in the solution of many 
kinds of boundary value problems [4,5,6]. In the construction of these circuits one con- 
sideration enters that is not present when the finite-difference equations are solved by 
purely numerical methods, namely that the circuit must be physically realizable. If the 
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circuit is to be constructed of resistors only, it must contain no ‘“negative”’ resistors, 
and the resistors must have the same resistance looked at from either end. In terms of 
the matrix of coefficients of the finite difference equations, the necessary and sufficient 
conditions that a network of resistors satisfying the equations is physically realizable 
are that the matrix must be symmetrical; that the nondiagonal terms in any row must 
all have a sign opposite to that of the diagonal term; and that the absolute value of their 
sum must be less than the absolute value of the diagonal term. These assertions can be 
easily verified by an examination of the equations of a network of resistors. Such re- 
strictions are not imposed on purely numerical solutions of the difference equations. 
Since the author of this paper is primarily interested in the solution of boundary value 
problems by means of electrical analogy, these restrictions have been imposed on the 
methods to be presented. These restrictions will have the effect of eliminating difference 
operators of higher than the second order in the equations for the network. 

The problem of changing cell size within a given rectangular gridwork has been solved 
by Southwell (Ref. 1, pp. 98-100) by a method which leads to network elements which 
are “physically unrealizable’”’ according to the rules laid down above. A network for the 


solution of Laplace’s equation 


re] “o 
Ox Oy” ° (1) 


is shown in Fig. 1. The finite difference equivalent of Laplace’s equation for non-ex- 
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Fic. 1. Network for changing cell-size in a network of squares by Southwell’s method. 


ceptional points in a square gridwork is, neglecting fourth order and higher difference 


operators, 
Lo Gn — 0) = 0, (2) 


where ¢, is the value of the function at the point in question, and ¢, is the value of the 
function at one of the nearest neighbors of the point in question. 

This equation is made to apply to exceptional points by placing upon it the following 
interpretation: the value of the wanted function at any point is equal to the average of 
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its values at the vertices of a square, the center of which is at the point in question. For 
example in the network of Fig. 1: 


dius + Po + oth — 45 = 0, 
dis + dis + 6 + os — 4610 = 0, (3) 


dois + bio + os + os — 4d, = (0. 


The network resulting from these equations is physically unrealizable because the 
matrix of coefficients is not symmetrical. For instance, the coefficient of ¢,. in the third 
equation is +1 but the coefficient of ¢, in the second equation is zero. 

The treatment that has been given to the problem of a curved boundary by writers 
on relaxation methods is largely intuitive (Ref. 1, pp. 67-78). From a consideration of 
the equilibrium of his tensioned network of strings Southwell concludes that, in the 
neighborhood of a boundary where the wanted function is known, the tension of a string 
connecting an interior point to the boundary should be inversely proportional to the 
length of the string. As applied to Eq. (2) this means that each term for which ¢, normally 
would fall outside the boundary should be weighted inversely as the length of the distance 
between ¢, and the boundary. 

The treatment of boundaries along which the normal derivative of the function is 
specified is less simple. In terms of the physical model, the transverse load at the edge 
is replaced by a statically equivalent set of forces applied to nodes just inside the boundary 
and to ‘‘ficticious” nodes just outside the boundary as in Fig. 2. The coefficients of the 
terms of Eq. (2) involving “‘ficticious’”’ nodes have values between zero and one. In some 
instances the coefficient for an element between “ficticious” nodes (e.g. the element 
between nodes 2 and 3 of Fig. 2) is set equal to zero and in other instances it is set equal 
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Fig. 2. Network of squares near a curved boundary. 


to 1/2. In a recent article, which was largely concerned with the development of accurate 
network formulas, the conclusion was reached that “‘the writer has failed to find any 
completely satisfactory method of dealing accurately with boundary conditions [in- 
volving a derivative] when the direction of the normal cannot be identified with that of 
a mesh line as in the case of curved boundaries”’ [7]. 

The difficulty with a curved boundary is caused by the fact that for a network of 
regular polygons the location of node points is unalterable. By using irregular polygons 
it will be possible to put the node points on the boundary. In this paper the situation 
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will be further clarified by giving a precise physical significance, in terms of the original 
field problem, to the terms of the resulting generalized difference equations. 

Derivation of the asymmetrical network. The problem at hand is the solution of the 
following equation together with appropriate conditions on ¢ and its normal derivative 
at boundary points of a finite plane region. 


V -(cVo) +7 = 0. (4) 


The following physical interpretation may be made of the symbols appearing in this 
equation: ¢ is the electrical potential in a plane region of conducting material. ¢ is the 
conductivity of this material. r is the density of currents inserted into the region from 
external sources. —oV¢? is the vector density of currents flowing in the material. ¢ and 
7 may be scalar functions of position and 7 may also be a linear function of @. 

This physical interpretation will aid in visualizing the constructions to be made. 

The problem of forming an asymmetrical network whose equations will replace 
Eq. (4) can be stated in the following manner. Given a region in which Eq. (4) holds and 
a large number of points in the region chosen at random, in what way should the points 
be interconnected with ‘physically realizable” electrical resistors in order that the 
voltages at the nodes shall be as nearly as possible the correct solutions of the boundary 
value problem characterized by Eq. (4) and appropriate boundary conditions? 

A unique answer cannot be given to this question at this time. A reasonable necessary 
condition that should be applied to the network is that for a homogeneous conductivity 
(¢ = constant) and a uniform field (V@ = constant), which is however arbitrarily 
oriented, the voltage at the nodes should give the exact solution of Eq. (4). It will be 
shown that more than one network connecting the given points can be constructed that 
satisfies this condition. 

In the method of solution that has been chosen the first step is to connect the randomly 
chosen points by a network of triangles, as in Fig. 3. The network should be planar (no 


Fic. 3. Asymmetrical network of triangles. 





cross-overs) and none of the interior angles of the triangles should be obtuse. It may be 
necessary to insert a few additional points in order to fulfill the last condition. 
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Consider a portion of this network shown in Fig. 4. The perpendicular bisectors of 
the sides of the triangles divide the region into polygons surrounding each point. A net- 
work of resistors is now constructed connecting the vertices of the triangles. The voltage 
across each resistor shall be interpreted as the line integral of the gradient of the potential 
between the two points it connects. For example, for points A and B of Fig. 4: 


(5) 








Fic. 4. Portion of the asymmetrical network of triangles. 


The current in the resistor shall be interpreted as the total normal flux crossing the 
common boundary of the dotted polygons surrounding the two points. Since the current 
density is —0V¢, we have: 


a2 


lan = —| o(Vo-n) ar, (6) 


1 


where n is a unit vector normal to dr. If V@ and ¢V@ are now expanded in Taylor’s 
series about the point 0, the midpoint of the segment A—B, and all terms except the 
first are neglected, Vé & (VW), and ¢V¢ & o.(V¢). . Since the segment AB is normal 
to the segment 1-2, the projection of V@ on dl is the same as the projection of V¢ on 


n. Hence 
Ve — Va = las | Vo |o cosa, (7) 
Tan ZS —Gri2 | Vo |p COS a. (8) 
The value of the resistor connecting A and B is 


V,-—V l 
AB Oo! 12 
Hence R,, depends only on the physical properties of the material and the manner 
in which the region is subdivided. If the segment 1-2 were not perpendicular to AB, the 
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value of the resistor would depend on the orientation of the field. It can also be shown by 
a simple geometrical argument that 


; 1 
ime” 


where 8, and £, are the interior angles of the triangles subtended by the segment AB. 
If both 8, and 8, are acute angles R,, will be physically realizable. 

In addition to calculating the value of R,, it is necessary to decide on an area element 
to be associated with the inhomogeneous term, 7, of Eq. (4). If Eq. (4) is integrated over 
the polygon 1-2-3-4-5 surrounding point B of Fig. 4, 


= * (etn 6, + etn p.), (10) 


I V7 -(cV¢) dS + I yim @. (11) 


B B 


By Gauss’ integral theorem: 


I V -(cV¢) dS = $ o(V¢-n) dr. (12) 

Hence . ‘ 
$ o(V¢-n) dr + I rdS = 0. (13) 

8 B 


If the surface integral in Eq. (13) is replaced by the value of r measured at B multiplied 
by the area of the dotted polygon and the line integral is replaced by network currents 
from Eq. (6): 

> | en + TaAp mend 0, (14) 


where J,, is the current flowing into node B from the p-th adjacent node. Eq. (14) is 
Kirchoff’s law for the sum of the currents entering a node. It shows that the appropriate 
area for calculating the current to be inserted into node B is the area interior to the 
dotted polygon surrounding B. 

By substituting Eq. (9) into Eq. (14) the generalized difference equation for node 
B is obtained: 


bs oo(782)(V, — Vz) + raAz = 0, (15) 
Pp Bp 
where /;, is the distance between node B and node p and rz, is the length of the segment 
that is common to the polygons surrounding node B and node p. 

The method that has been described will work for any network configuration in 
which the perpendicular bisectors of the branches meet at a point. Thus besides for 
triangles, the method will work for rectangles, regular hexagons and isosceles trapezoids. 
Later it will be shown that the perpendiculars to the sides need not bisect the sides, so 
long as they meet at common points. 

From the method of derivation given here nothing can be inferred as to the accuracy 
of the solutions of Eq. (15), except that for a region of uniform conductivity with a 
constant, arbitrarily oriented potential gradient, the solutions will yield correct answers 
to the field problem. The magnitude of the errors will be investigated in a later section. 
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Applications of the asymmetrical network. The manner in which the asymmetrical 
network can be applied to the problem of a curved boundary is illustrated in Fig. 5. 
A certain number of points are placed on the boundary and lines joining boundary points 
are considered in the same manner as lines joining interior points, except that the con- 
ductivity of material outside the boundary is set equal to zero. If the outward normal 
gradient of the field, 6/dn, is specified at the boundary, an additional current equal to 
d¢/dn multiplied by the conductivity and the length of boundary associated with each 





‘ 


Fic. 5. Asymetrical network near a curved boundary. 


node is fed into each boundary node. For point A of Fig. 5 for example this additional 
current is o4[8¢/dn], 712. This current is equal to the total flux crossing the boundary 
along the segment 7,2 . Hence for the asymmetrical network, boundary points are treated 
in almost the same manner as interior points. 

The manner in which the principles of the asymmetrical network can be used to 
change cell size in a network of squares is illustrated in Figs. 6a and 6b. The numbers 
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Fic. 6. Two ways in which to double cell size in a network of squares. 
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beside each branch of the network are the ratios r;,/l,, from Eq. (15). For example the 
branch connecting nodes B and C of Fig. 6a has a length a/ V2 while the length of the 
common boundary is a/2‘V 2 giving to rgc/Igc the value 1/2. It will be noted that the 
only values of this ratio occurring in these figures are 1/2 and 1. In contrast with the net- 
work of Fig. 1, the networks of Fig. 6 can be physically constructed and used for an 
analog computer solution. Another advantage that should not be underestimated is 
that the current in each element of the networks has a real physical significance. It 
represents the total normal flux crossing a known line segment. 

An example of the application of the asymmetrical network to a complete problem 
is illustrated in Fig. 7. The problem concerns the calculation of the resonant frequencies 








and field patterns of the so-called conical-line cavity resonator, a cross-section of which 
is illustrated in Fig. 7c. For TEM modes (modes in which electric field lines lie in planes 
passing through the vertical axis and magnetic field lines are concentric circles surround- 
ing the axis) the equation governing the variation of the magnetic field in a plane passing 


through the vertical axis is 
? ie 
V-\- VAH;) + — As = 0, (16) 
p p 


where p is the perpendicular distance to the vertical axis, H; is the covarient component 
of magnetic field intensity and )” is an eigenvalue related to the frequency of oscillation. 
The physical component H, is equal to H;/p. The boundary condition on H, is that 
0H,/dn = 0 along the walls of the cavity. The cavity is assumed to have the shape of a 
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ANALYTICAL SOLUTION 





Ha-(40.0 





Fic. 7b. Computer solution of the lowest mode of a conical line resonator. Measured wavelength = 
3.984a. Correct wavelength = 4a. 


sphere where 30° conical dimples at the poles. By comparison with Eq. (4) we see that 
o = 1/p andr = )*H;/p. 

Solutions for this problem were obtained by the electrical analogy method. An 
electrical circuit together with numerical values of the network elements are shown in 
Fig 7a. Since the coefficients of the second term in Eq. (16) is inherently positive and 
varies with the frequency of oscillation, a variable “‘negative’’ resistance is required for 
its realization. This difficulty is avoided by using inductors for the elements between 
nodes and capacitors for the negative elements. If the network is resonating with a fre- 
quency w, the equation for the sum of currents at any node, B, is: 


I 


p wl gy 


(V, — Vz) — twlgV,z = 0. (17) 


Upon multiplying this equation by it and comparing the result with equation (15), 
using the values of ¢ and 7 appropriate to this problem, it becomes apparent that: 


Ls, = p ize Cp =pAp, w =D. (18) 
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57a. 


1.4 


TEM mode of a conical line resonator. Contours are lines of constant H3, which are parallel to the electrical 
field. Measured wavelength 


r 


Higher 





7¢. 





Fia. 





The resonant frequencies and the corresponding eigenvalues were obtained experi- 


mentally. 
The eigenfunction for the lowest mode of the conical-line resonator has the simple 


form 


H, = sin —, (19) 
2a 
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where r is the distance from the center and a is the radius of the sphere. This solution is 
compared with values measured on the network in Fig. 7b. Lines of constant H; , which 
are parallel to the electric field, are plotted in Fig. 7c for a higher TEM mode. In Fig. 
7a it is seen that the principles of the asymmetrical network have been used to fit the 
location of nodes to the natural boundaries of the cavity resonator and also to effect 
changes in cell size so as to keep the cell area nearly constant in going from the center to 
the outer wall. The solutions were obtained on the California Institute of Technology 
Electric Analog Computer [8]. 

The use of Taylor’s series expansions. A network for changing cell size within a 
network of squares, which has not been constructed according to the principles of the 
preceding section, is shown in Fig. 8. This network appeared without much explanation 
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Fia. 8. Method of Reference (9) for doubling cell size in a network of squares. 


in a paper by Spangenberg, Walters, and Schott [9]. The numbers beside each branch 
are coefficients corresponding to rz,/ls, in Eq. (15) and the numbers beside each node 
are the areas to be associated with each node. If this network had been constructed by 
the methods of the preceding section there would be, for example, a branch connecting 
nodes A and B. The values of the branch coefficients shown in the figure can be obtained 
if perpendiculars are drawn not through the midpoints of the branches but in the manner 
shown in Fig. 8. (These lines were not drawn in the paper quoted above.) This is per- 
missible since, in the derivation leading to Eq. (15), the assumption that the perpen- 
diculars to the branch segments were drawn through the midpoints of these segments 
was not essential. If the area of the dotted polygons surrounding each node are now 
computed, the results do not agree with the values given in the figure. For example the 
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area of the polygon surrounding node C is 3/4 a° rather than 5/8 a’. The question is thus 
raised as to which of these values is better. This and other questions can be investigated 
by means of an error analysis using Taylor’s series expansions. 

The following investigation is limited to the case where the coefficients of Eq. (4) 
are constants, and for simplicity these coefficients will be set equal to unity. The results 
of the analysis are valid for any values of r and o. With 


V*a+1=0, (20) 


the corresponding generalized difference equation for the potential ¢) of any node of 


the network is 


> Yi, — $0) + Ao = 0, (21) 


where @, is the potential at a neighboring node and A, is an element of area to be associated 


2 


$y a comparison of these equations we see that a measure of the error introduced 


with @ . I 
by replacing Eq. (20) by Eq. (21) is 

é& = Ao(Vb)o — D> Yolo — o)- (22) 
If the point where ¢, is defined is taken as the origin of a cartesian system of coordinates, 
the potential at any neighboring point can be obtained by a Taylor’s Series expansion. 


The summation in Eq. (22) can be expressed in terms of such a series as 


} > l - 0 0 \" 
P — = < 4 + y 
. V wold $0] =n al ” Ox Yp 2) ?. (23) 


0 


The coefficients of the terms of this series are arranged in tabular form below: 


term coefficient 

/ fel 1 . > 

| 4y ) ( i » 2 } pot p 
Od v 


( 46) C, = r& Y oY 


\ay/ > 
Ip a 
- C; = Y ox? (24) 
\ 94°" 9 pop 

0d ‘ 
(. 1, ) ( Sees z2 } vot pYp 
\or OY ‘ Pp 

0 aries 
(72) C.=4D You? 
\dy"/ 5 = 

ete. etc. 


In this table all coefficients except C; and C; should vanish independently if the 
finite difference approximation is to be correct. The vanishing of the first two coefficients 
is equivalent to the statement, pertaining to an analogous problem in statics, that the 
center of gravity of loads, Y,. , concentrated at the surrounding node points should be 
at the origin. It will be demonstrated that this will be the case if the Y,o’s are calculated 














1953] AN ASYMMETRICAL FINITE DIFFERENCE NETWORK 307 


by a simple geometrical method to be described. In Fig. 9 a closed polygon is shown 
whose sides are each drawn perpendicular to line segments radiating from a common 
point. A concentrated load is placed at the end of each line segment equal to r,/l, where 
l, is the length of the segment and r, is the length of the side of the polygon perpendicular 


Ye Y.* = 


rd 
i t 
a“ 





Fic. 9. Construction to prove that coefficients C,; and C, vanish. 


tol, . The moment of the loads about any line, X — X, drawn through the common point 
making an angle a, with each radiating line is 


Tp 


Moment = pe 


Pp 


l,sne, = >. 7,sina, (25) 
P p 
The last expression is the sum of the projections of the sides of the polygon on the line 
X — X, which vanishes because the polygon is closed. Note that in this proof the polygon 
need not be drawn through the midpoints of the radiating lines. 

In general the coefficients C,; and C; of Eq. (24) will not be equal. If they are not, 
these terms should be combined to give a Laplacian and a “hyperbolic” operator, i.e., 


F *s) (2$) C, +. f, (zs ¢) c da C; ($ 7) 
ci—> Ci—s) = 2 (5-4 4+ 3 —3__—-8 (-* _ —4). 
(2s 0 sites dy / > 2 Ox” - dy’, + 2 Ox” dy"/ > (26) 
The second term on the right can be put into the form, (0°¢/d£ 9). , by a coordinate 
rotation through 45° while the first term is invariant in a coordinate rotation. From 
Eq. (22), the coefficient, (C; + C;)/2, should be equal to A, , and this provides a means 
for choosing the appropriate area for each cell. 
C; it 1 —< > 
Ay = SE“ = 2D Vela + y). (27) 
2 i:‘sS 
If Y,o is again considered as a load acting at x, , y, the appropriate area should be 1/4 
of the sum of the polar moments of inertia of the loads. 
If the value of Y,. is substituted into Eq. (27), then, 


1 
Ao = 4 p> r,l, 
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since 
2 2 2 
L, = I, T Yo « 


If the dotted polygon of Fig. 9 passes through the midpoints of the radiating line 
segments, the area of the polygon satisfies Eq. (28). In the general case it seems desirable 
to construct boundaries which clearly delineate the area to be associated with each node. 
A construction which satisfies this condition and Eq. (28) is shown in Fig. 10. The ap- 








Fic. 10 Construction for obtaining the region appropriate to each nodal point. 


propriate area is the area enclosed by joining the midpoints of the radiating line segments 
to the points of intersection of the perpendiculars to the segments. This procedure gives 
the same values of area coefficients for the network of Fig. 8 as those computed by 
by Spangenberg, Walters, and Schott [9]. 

Referring to Fig. 4, it can be seen that a network capable of solving the potential 
problem can be constructed in which the branches are the lines drawn perpendicular 
to the branches of the original network and in which the nodes are the vertices of the 
“dotted” polygons. For this new network the roles of r, and 1, are interchanged, and 
this new network may be called the “geometrical dual’’ of the original one. Any network 
constructed according to the rules given in this paper has a geometrical dual. It will be 
observed that Fig. 8 is very nearly the geometrical dual of Fig. 6b. 

The values of Y,. and the area coefficient A, have been chosen so as to satisfy three 
conditions at every node. Since on the average there will be about five or six branches 
radiating from each node it seems that it should be possible to choose the values of the 
Y,0's to satisfy about three more conditions. However each Y,o occurs in the equations 
for two nodes so that on the average for each node there are only about three independent 
coefficients. The non-vanishing terms of Eq. (23) will be simply regarded as errors which 
can, if desired, be incorporated as correction terrfis in the late stage of the calculation. 
€) in Eq. (22) can be estimated from a trial solution and added to A, in Eq. (21). This 
follows the general procedure used in relaxation calculations [3]. 

The error series can be easily calculated from Eqs. (22) and (23) for each node of any 
given configuration. For a network of squares, with branch length a, the terms of lowest 
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order in the error series are 


_ _a (a 2¢) 
— ae (7s dy*/,° (29) 


For a network of equilateral triangles with branch length a, 
a‘ 
= 16 (V*d)o ° (30) 


For nodes without double symmetry in their branch patterns terms of lower order will 
occur. 

The following table gives the errors at each type of asymmetrical node in Figs. 
6a, 6b and 8. Terms of the fourth order and higher are not included, 














Fig. 6a " 

Nate A Sly -aa)t 

Node B E (4 - =) ? 3 (75) 12 

Node C | -< (4% ar 2) ? (595) | 

Fig. 6b , 
cues (£(8-)-RGFs 28] 
Node B | -% (3 7 *) + + (sa) | 

Fig. 8 

Node A E 5 + a (3 + Er = ay Ai =) | 
Node B | -¢ 5 + (-5 . say - ar a 2)\ 9 
Node C [$5 (astay - Fay) | 


For each figure, the sum of the e)’s in a vertical strip of width a is equal to (a*/8)(d*¢/dy’*). 
This iridicates that the other terms contribute to a merely local distortion of the field 
pattern and that on a large network this distortion will be negligible. It is interesting 
that for a one dimensional network in which a two to one change in cell size is made, 
the leading term in the error series at the point where the change is made is just (a*/8) 
(d°¢/dy’). 

Conclusion. A method has been described for constructing an asymmetrical finite 
difference network that can be used in the solution of second order boundary value 
problems. The coefficients of the difference equations that govern the network can be 
found by simple geometrical measurements. The asymmetrical network has the 
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advantages that it provides a simple solution to the problem of fitting a gridwork to a 
curved boundary, and that it provides a means of changing cell size in such a way that 
the network is “realizable’”’ by means of physical electrical elements. A clear interpreta- 
tion has been given to the currents which flow along the branches of the network. 


. R. H. MacNeal, The solutior 
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LEGENDRE FUNCTIONS OF FRACTIONAL ORDER* 


BY 
MARION C. GRAY 
Bell Telephone Laboratories, Inc., Murray Hill, N. J. 


Introduction. In the theory of the propagation of spherical waves in free space 
the angular wave functions are Legendre polynomials, P,,(cos @), or associated Legendre 
polynomials, P", (cos 6), where n and m are restricted to integral values. These func- 
tions are polynomials in cos 6, their properties have been widely studied, numerical 
values have been tabulated, and in general they may be regarded as known functions. 

In more recent years, however, Legendre functions of non-integral order, which we 
shall denote by P,(cos @), have also occurred in physical problems. Thus, for wave 
propagation inside a circular horn of given angle, the boundary conditions introduce 
a characteristic equation which is actually an equation in the parameter v. It has been 
customary to simplify the problem by choosing horn angles corresponding to integral 
values of vy, but a complete solution should include a study of the behavior of P,(cos 6) 
as a function of v. 

Similarly, in the mode theory of antennas developed by Schelkunoff the appropriate 
angular wave functions in the antenna region are Legendre functions of order n + 120/K, 
where n is an integer and K is the characteristic impedance of the biconical antenna 
to the principal wave. For thin cones K is large and the order of the Legendre functions 
is nearly, but not quite, integral. Further, when the cone angle is large, »y may have 
quite general real values. 

Another application has appeared early this year, when P. Grivett used Legendre 
functions of fractional order in the approximate solution of an electron lens problem, 
with particular emphasis on small values of »v. 

Thus it appears that the properties of Legendre functions of non-integral order 
are of quite general interest, and it may be worth while to put on record some formulas 
that were developed a few years ago in connection with Schelkunoff’s antenna theory. 
At that time the formulas were used to compute values of P,(cos 6), 0 < 6 < a, for 
values of v between 0 and 2 at intervals of 0.1, and curves based on these computations 
have already been published.** Those curves show P,(cos 6) as a function of 6 for the 
fractional values of v; in this memorandum we include a table of numerical values (Appen- 
dix, Table I), and also a new set of curves (Figure 1) showing P,(cos 6) as a function 
of v for values of @ between 0° and 175°. We have confined our computations to real 
values of v, but it might be worth noting that the approximate formulas, and in particular 
the fundamental series expansions (3) and (17), are also valid for complex values of v, 
in all regions in which they converge. 

The function P,(cos @) has a logarithmic singularity at @ = 7 for all non-integral 
values of v; and it may be expressed in closed form at 6 = 2/2 for all values of v; hence 


*Received September 8, 1952. 

tP. Grivet and M. Bernard, Théorie de la lentille électrostatique constituée par deux cylindres coariaur, 
Ann. de Radioélect., 6, 1-9 (1952); P. Grivet, Un nouveau modéle mathématique de lentille électronique, 
Jour. de Phys. et le Rad., 13, 1A-9A (1952). 

**S. A. Schelkunoff, Applied mathematics for engineers and scientists, D. Van Nostrand, New York, 
1948, pp. 423-424. 
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P, (cos 0) 





Fia. 1. 


it is convenient to consider separately appropriate expansions in the neighborhoods of 
6 = 0, x/2, x respectively. 

It may also be pointed out that for nearly integral values of v, say v = n + 4, it is 
sufficient to consider the case of n = 0 and small (positive or negative) values of 6. 


Then the recurrence formulas 


1 + 26 6 , : 
P,.3;(cos 6) = —— cos OP;;(cos 6) — ———~ P_,(cos 8), 
ve Ss eS ) Os @) i+ s (cos 6) 
(1) 
2n — 1+ 26 ’ a-i+s, ; 
P.,;(cos 0) = — : + _— cos 6P,,-;.;(cos 6) — — : + . os P,,-2+3(cos 6), 
may be used to obtain the required values for all values of n. 
1. Small values of 6. When @ is not too large, the usual series expansion,* 
, —(—)‘(v+s)! .,, 1 
D Sain nie ee tee iG) 
P,(cos 0) = >> > = aed sin” 5 6, (2) 


2=0 


converges quite rapidly and tables of the factorial function are available. This series, 
however, does not exhibit the analytic nature of P,(cos @) as a function of v. Hence, 





*S. A. Schelkunoff, loc. cit., p. 420. 
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for small values of v, we developed a series expansion in powers of », 


P(cos =) = 1+ > a,’, (3) 
n=1 
are functions of @ which can be computed from a set of recurrence relations. 


where the a’s 
= sin’ (6/2) we can express the a’s in the following form: 


If we write z 


) ecie oo 2" (-)**? i) 2 
Qonsi = r k..- Lang =O k... 3 
— n! p> ake cna a =<, "Ss" (4) 


where 


(5) 


ee s=n+1,n+2,-:- 


It can be seen that the values k can be tabulated very rapidly, and then multiplied by 
the appropriate factors involving z to obtain the values a. 


In particular 


a, = ->- = log (1 — 2) = 2 log cos 5 6, 
s=1 * 

i (6) 
s=1 8 


8 
~ 
~ 


x ~* ro) e-l 1 
~ 
a = Zz 2.0 a= oe G2, 2) a." p> 2° 
s=2 s ‘ s=2 r=1 7 


s 


When @ S 2/2 we have z S 1/2 and the series (4) converge quite rapidly, while the 
successive a’s become smaller. The series (3) is valid for either positive or negative 
values of v, and can be used very conveniently to compute P,(cos @) for values of »v 
such that | vy | < 1/2. Then the recurrence relations (1) may be used for larger values 
of v, using also the general relation P_,_, = P,. 

At vy = +1/2 the convergence is rather slow, except for small values of 6, but the 


Legendre functions can be expressed in terms of elliptic integrals, 
oe . 

P,,.(cos 6) = 2E\sin = 6} — K\sin = @] |, 
7 2 2 


P_,/.(cos 0) = ? K(sin } a), 


(7) 


where K and E are the complete elliptic integrals of the first and second kind, respectively. 
These functions have been frequently tabulated, and may be regarded as known. The 
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recurrence formulas enable us to express the Legendre functions of order »y = n + 1/2 
in terms of K and E£. 
We shall see later that, for small values of v, the first approximation 


l 
P,(cos 8) = 1 + 2» log cos 5 0 (8) 


is good throughout the range 0 S 6 < r. 
2. The neighborhood of @= x/2. When 6 = 7/2 the value of the Legendre function 


I cos 4vm(4v — 3)! 
P,\ cos 57) = (9) 
“ V r (3)! 


for all values of v. In the neighborhood of 7/2 we write 6 = 7/2 — a and a series expan- 
sion which converges rapidly for small values of @ is 


1 > 1 Reo: acs, BYE ae 
l : (sv + 47 x)! v r (2sina) 
Pj cos e r= a) = P,(sina) = » " ae , COS 5 : 


. 
0 Var (fv — 3r)! “ rs 


1 (}» — 4)! ,/1 l = er 
= COS 5 ur =— p y+ > 9% 95m a 


- VU r (3)! \2 - ~ 
, 9 (4y)! {1 l - ee 
+ 2sina sin 5 vr — = — + oe ~se oa 535) SIN a}. (10) 
2 Var (fv — 3)! \2 2 2’ 2 


When we consider P,(cos 7/2) as a function of v, it ean be shown that the first few 
terms in the expansion of the function (9) in powers of v are 


-\ e rit. 2 
P,( cos 3) = 1—vlog2 — 1 E nr — (log 2) | 


v ul eo? — 93 _ 3 5 ] Gite 
+2] tog2 (log 2) = de + (11) 


ae | 
n=} Mt 


These terms. however, check with those obtained in Section 1, since, at 0 = 7/2, 


l l 
2 log C08 5 6 = 2 log cos= ax = —log 2; 


4 
the relation 


= ] r + : 
~. a © — = (ing 2 
2 5 = 12 ~ 2 “08 


is given in Smithsonian Mathematical Tables, p..142; and 


il 


“02,5 “i p l : 
2 = = — (log 2) — 9 (log 2)° — 4 27 


2's 12 
can be verified numerically. If we try to expand the series in equation (10) in powers 
of v we find that when all terms have been collected we simply return to the series (3). 

3. The neighborhood of 6 = =x. In this region the most useful expansion for all 
values of v seems to be that first obtained by E. Hille.* If we write 6 = r — ¢ Hille’s 


*See E. W. Hobson, Spherical and ellipsoidal harmonics, University Press, Cambridge, 1931, p. 225. 
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formula may be written in the alternative forms 
P, [cos (x — ¢g)] = P,(—cos ¢) 


sin v — (—)'(pv ‘)! . A 
_ sin vm > | "( +o g log sin 5 @ + Yo ++ u(—rtr— 0-240 | (12) 


T 0 (vy —r)! 


2 sin yr — 
= | log sin 5 ¢ + cos vm |P,(cos ¢) 


us “ 


r 


, sinvr = (—)'(v +7)! oes 
fe cpr We tN + He -— 9 - WT 13) 


where z = sin’ ¢/2 and (x) is the logarithmic derivative of the factorial, as used by 
Hobson, 


d 
x = °! 
¥(x) de (log x!). 


When ¢ is small the series in equation (13) converges rapidly and only a few terms 
are significant. Further, the first term of the series vanishes with v for small vy and thus 
the dominant terms in the expansion in powers of v are those obtained from the first 


term in equation (13), 
— 
P,(cos 0) = P,(—cos ¢) = P,(cos ¢g)\ 1 + 2» log sin 5 ¢ 
(14) 


~ 1+ 2p log cos 5 6. 


If we include the series terms, the expansion in powers of v for small g found most con- 


venient for computation may be written 


P,(—cos ¢) = P,(cos a cos ve + 25m or (Iog sin ; eg + Wr) — vo) | 


sin vr 2 - 
a (Co + cw + cy” + --°), (15) 
where 
fc = Gs, 
Cc, = —2a,, 
2’ 
Co = —a,; +2 2d, 
ii (16) 
@ 2° 
Cc = — 2a, — 2 > o2.. —; 
a=2 s 
= 2° ~ 2" 
¢; = —a, —2 Do. 3 — 2 Do 93, 3 
a=2 8 a=2 § 
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and where the values a are the constants of P,(cos g), z = sin’ g/2, and 


s—l 1 


On,s = ne 
rai I 


When ¢ is not too small it is also possible to use the series expansion 


P(—cosg) = 1+ >> by’, (17) 
n=1 
with 
= 
b, = 2 log sin 5%) 
r? 
b, = a,b, —a—- > 
Vo a0; a2 6” 
(18) 
2 a s es) 1 
b, = b(a, - =) +2>5-2D5, 
0 a=] § s=1 s 
© — | r r* = 2° 
b= bi(ay — Sa) — 2m D+ Ea — a+ Be - 2 Dawe, 


The coefficients are of course more complicated than those of equation (3), but if the 
a’s have already been computed the remaining terms may be evaluated without too 


much labor. 
4. Zeros of P,(cos @). In Grivet’s electron lens theory certain focal distances are 


determined from the roots @ of P,(cos 6) = 0 and from the values of dP,/dé@ at @ = @. 
From the values of Table I the curve of Figure 2 has been drawn, showing the roots 


“a Ne ; an 


80° 
ZEROS OF R, (cos 6) ze 
© APPROXIMATE VALUES ae 
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Fia. 2. 


of P, for values of v between 0 and 2. However, approximate values of the roots can 
be found from the approximate formulas of the preceding sections, and it is surprising 
that the simple formula (8) gives values very near the true values except in the neigh- 
borhood of » = 1/2. At this point the root can be obtained from elliptic integral tables. 


For smaller values of v we find from (8) 


1 1 , 
C08 5 6 = exp (-1), (19) 
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and for values of »y between 0.5 and 1.5 we combine equation (8) with the recurrence 
relation (1) to derive the equation 
1 1 (1 + 26) cos 6 — 6 
les cn = (2 -— yv=1+ 6. 20 
B89 25 (1 + 28) cos 6 + 8’ + (20) 
Similar equations can be obtained for larger values of v by repeated use of the recurrence 
formula, and it is always possible to obtain an equation for the root which expresses 
log cos 6/2 as the ratio of two polynomials in cos 8. 
When » is small the root 4) is near x, and a somewhat more accurate formula is 
obtained from equation (15): 


= oie log cos 49 T 
lor sin = @ 02 mee on = ~ 

ogsin5¢ = Ty oy iconts 3 cot »r — ¥(v) + ¥(0), (21) 
where 6 = x — ¢. Similarly for »v = 1 + 6, where 6 is small, the smallest root can be 
found from equation (20), but there is a second root near + which is determined more 
accurately from the equation 


(1 + 28) cos 6 — 5] _ cos 6[y(8) — ¥(0) — log sin 36] (22) 


(1 + 28) cos 0+ 6 (1 + 28) cos 0+ 5 


In Figure 2 we have indicated by circles the values of the roots obtained from equa- 
tions (19) and (20) where these may be distinguished from the curve values. 

For the derivative of P, at @ = 6) we can find simple formulas by differentiating 
the approximate formulas (3) and (15). Thus retaining the first three terms in (3) and 
using the approximation (19) for 0 we find 








] 
log cos = 


T 
_- 5 cot in| 


dP, | 2v (23) 


dO \pnn, SiN O 





When » is small, so that 6 is near x the asymtotic approximation is 


aP, 2 sin yr 1 (24) 





d6 |eno, | 7Sin 0) 1 + 2 logsin 6/2 
At v = 0.5 the derivatives of the elliptic integrals give 


da | — —K(sin 60/2) 
6 P,,.(cos 6) , von & (25) 


as shown by Grivet. For v > 0.5 we may combine (20) and (3) with the recurrence 
formula to find approximately 





d _ = 28(1 + 8) (26) 
dé 0=85 sin 6[(1 + 25) cos 8 ++ 5] 
Note that this equation is valid as 6 — 0. For we have P,,; — cos 6, and the limit is 
approached in such a way that cos 0) = lims,.. 5/(1 + 26). Thus in equation (26) the 
limiting value is 


P,,;(cos 6) 





| 
1 
0=86 sin 0 





dy 
doe i(cos @) 





which is correct. 
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110° 880955 
120° 852374 
130° 817704 
140° 774511 
150° .718190 
160° 638358 
170° .501717 
175° 365201 
] y =1.1 
10° .982463 
20° 930510 
30 846085 


40° 732358 


50° 593260 
60° 435048 
70° 262707 
80° +.033160 
90° —.096662 
100° — .269685 
110° —.423902 
120° —.567487 
130° —.679240 
140° —.756377 
150° —.792525 
160° —.775459 
170° —.672324 
173° —.542651 





> 


998171 

992665 
983428 
970362 
953322 
932102 
906416 
875872 
839927 
797813 
748422 
690081 
620144 
534092 
423320 
268268 
005894 


254581 


979971 
920782 
825099 
697230 
542960 
369302 
184212 
003726 
185629 
352628 
496128 
607957 
680298 
705113 
672251 
562976 


318188 


058480 





Le ge ndreé 


3 


997028 
988095 
973140 
952059 
924694 
890814 
850092 
802069 
746089 
681210 
606031 
518406 

414869 
289416 
130467 
088558 
453932 


813813 


977330 
910521 
803106 
660787 
491146 
303277 
107352 
085870 
265506 
421019 
542575 
621232 
648801 - 
617004 
514788 
319350 
047601 


413527 


.974544 


MARION C. GRAY 


APPENDIX, TABLE I 


Functions of Fractional Order, F. 


4 


995735 
982927 
961544 
931521 
892752 
845072 
788227 
721889 
645288 
557670 
456324 
342882 
209624 
051166 
145688 
411669 
847492 
272544 


1.4 


899731 
780153 
623160 
438457 
237487 
032911 
162209 
335101 
473792 
567190 
608117 
587376 
498494 
331564 
063159 
395130 
833237 


—1 
—1 


— 


1 


.o 


994289 
977168 
948665 
908821 
857676 
795249 
721505 
636309 
539353 
430035 
307261 
169084 
012012 
170483 
391682 
683193 
150000 


599553 


971609 
888427 
756288 
584483 
385144 
172439 
038356 
231828 
393447 
510309 
571540 
570350 
500235 
337352 
135126 
186880 
696918 
165932 


—1 
-1 


T 


1 


6 


(co 
4 


992693 
970822 
934528 
884042 
819665 
741748 
650659 
546730 
430189 
301038 
157754 
002434 
170817 
366364 
596024 
888957 
343918 


774742 


1.6 


968530 
876616 
731556 
544892 
331534 
108625 
105747 
293895 
439820 
530315 
555254 
510520 
392577 
202235 
061645 
.412909 
930007 
383435 


990947 
963895 
919163 
857275 
778935 
685007 
576469 
454374 
319752 
173516 
016269 
151995 

332468 
528746 
749809 
019355 
420160 
791031 


-1 
—1 
—1 


965306 
864311 
706009 
504529 
277906 
046528 
168608 
347694 
473745 
533992 
520788 
431952 
270327 
042220 
+ 246321 

599730 
1.077683 
1.476543 


Ss 6) 


8 
989050 
956393 
902601 
828616 
735715 
625480 
499745 
360536 
209982 
050203 

116849 
289652 
467557 
651698 
847187 
069887 
378607 
652930 


-1 
—1 
—1 


1.8 


961937 
851520 
679698 
463536 
224546 
013395 
226363 
392684 
495011 
521924 
468913 
338567 
139860 
113682 
407651 
735204 
1.130627 
1.433916 


+. 


> 
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9 
987003 
948322 
884877 
798169 
692387 
563646 
421314 
266528 
102787 
— 066312 
237227 
406673 
574395 
734616 
885632 
041354 
227945 


378654 


—1 
—1 
-1 


1.9 


958424 

838256 

65267 3 

422057 
+ .171992 
070707 
278485 
428461 
503662 
495078 
402560 
234713 
007413 
260689 
536224 
8113058 
1.087486 


1 265476 


1.0 


984808 
939693 
866025 
766044 
642788 
5 
342020 
173648 
0 
173648 
342020 
642788 
766044 
866025 
939693 
984808 


996444 


2.0 


954769 
824533 
625 
380236 
119764 
25 
324533 
454769 


454769 
324533 
125 

119764 
380236 


.625 


824533 
954769 


989351 





THE LOCATION OF THE ROOTS OF POLYNOMIAL EQUATIONS BY THE 
REPEATED EVALUATION OF LINEAR FORMS* 


BY 
L. TASN Y-TSCHIASSNY 


University of Sydney, Australia 


1. Introduction. The author was recently engaged in problems connected with the 
solution of polynomial equations with the aid of an electrolytic tank analog’. In con- 
nection with this work he evolved a simple and apparently novel computational system 
of locating the complex roots of polynomial equations. This system is particularly 
suitable for ‘punched card’”’ and ‘digital electronic” computing machines, because 
it is essentially the evaluation of linear forms, repeated systematically. The present 
paper describes the principles of the system. 

2. The connection between a polynomial and a two-dimensional field. Let the 
polynomial the zeros of which are to be located’ be 


GZ) = S LZ = Sur + alyz (1) 


Let u be a scale factor and Z’ an auxiliary variable 
Z=uZ’ (2) 
As Lucas® pointed out, a rational function H(Z’) which has m > n first order poles at 


arbitrarily selected points Z,(s = 1, 2, --- , m) and whose zeros coincide with those of 
G(uZ’) can be derived by dividing kG(uZ’) by 





F(Z’) = [| (2 -— Z), (3) 
s=l 
where k is a constant. H(Z’) can be expressed as the sum of partial fractions 
ry = - oe) 2 F_A_ 
H(2') = k pon = p> 7-7" (4) 
where 
ee , 
A,=a, +2," =k BAZ’) (5) 
and 


a=(a+1) tom 


B(Z) = JI (@-2Z). (6) 

qa@=lto(s-1) 
As is the residue of H(Z’) at Z!. By eliminating in (1) to (4) the expressions Z’, G(Z), 
F(Z’), and H(Z’), we obtain two power series in Z. The comparison of the first coefficients 


*Received Nov. 7, 1952. 
1L. Tasny-Tschiassny and A. G. Doe, The solution of polynomial equations with the aid of the electro- 
lytic tank, Aust. J. Sci. Research 4, 231-257 (1951). 
*Capital letters stand for complex numbers, small letters for real ones. 
3F, Lucas, Résolution des équations par |’électricité, C. R. Ac. Sci. Paris 106, 645, 1072 (1888). 




























320 L. TASNY-TSCHIASSNY [Vol. XI, No. 3 


shows that 
>> A, = 0, if m>n-+ 2, 
(7) 


ku"""L, , if m=n+1. 


1 


The two-dimensional vector corresponding to the conjugate of H(Z’) is proportional 
to the field strength at Z’ of a two-dimensional field produced by sources of the intensity 
as” and vortexes of the intensity aS" positioned at the poles Z/ . If b Myer A, = 0, the 
field is self-contained, if b Bape A, ~ 0, a source and a vortex of intensities given by 
(—o1 A,) is to be imagined at Z’ = o. Since the zeros of H(Z’) and G(uZ’) coincide 
and the zeros of H(Z’) are saddle points of the potential, these saddle points determine 
the roots of the equation G(uZ’) = 0. This relation has been utilized for electrolytic 
tank analogs*”’. 

The second statement of (4) can be employed for a purely computational explora- 
tion of the field conditions. By a systematic cut and try method one can approach the 


points Z/ for which H(Z’) = 0, with any desired accuracy. In general, this method 
will be very cumbersome, the main reason being that the accuracy required in the 
computation of the terms A,/(Z’ — Z/) is considerably greater than the accuracy 


obtained in the result H(Z’). 
In the special arrangements discussed in the present paper the described computa- 
tional exploration becomes very simple, because it is essentially the evaluation of linear 


f= Dab, (8) 


in which the set of quantities a, depends on the special numerical problem in hand, 
and the set of quantities b, is taken from tables compiled once and for all. Digital elec- 
tronic and punched card machines are particularly suitable for this type of work, but 
it appears that a satisfactory efficiency can be obtained with ordinary commercial 
multiplying machines. When evaluating linear forms on these machines, no need arises 
to make a record of intermediate products, because when a certain number has appeared 
in the result register, a further multiplication adds or subtracts the additional product 
to this number. If the numbers a, are recorded on a strip of paper in a way that in a 
certain position all numbers a, can be made adjacent to the corresponding numbers b, 
of the table, the corresponding multiplicands and multipliers can be directly read off 
without risk of errors. The tables for b, can be arranged in a way that the same a,-strip 
can be used for different b,-sets. 

3. The computation of the residues A, for symmetrically arranged poles Z; . A 
convenient number is chosen for the scale factor u in a way that at least some of the 
points Z/ for which G(uZ?) = 0, are within the circle of unit radius with the origin 
as centre (in the following called the “unit circle’). Let the poles Z/ be the complex 
roots, of the equation 


forms 


Zz = 1, 





4A. R. Boothroyd, E. C. Cherry, and R. Makar, An electrolytic tank for the measurement of steady-state 
response, transient response, and allied properties of networks, Proc. Instn. Elec. Engrs. 96, 163 (1949). 

5A. Bloch, Solution of algebraic equations by means of an electrolytic tank, VIIth Internat. Congr. 
Appl. Mech., London, Paper No. IV-28. 
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1.€., 
Rafe, (s = 1,2, ---, m). (9) 
The quantity B,(Z!) [equation (6)] may be expressed as 
tua. & tl ah 
Bit) = tim | Fz ” 


which after substitution from (9) and evaluation of the indeterminate form leads to 
BAZ!) —_ me 27**/™ (11) 


From (1), (2), and (5) we obtain for k = m 


a=n 
A,= > (Lue, (@ = 1,2, >>: , m). (12) 
q=0 
The real and imaginary parts of A, are linear forms of the type of equation [8], because 
the trigonometrical functions appearing in them can be tabulated once and for all, if 
m is fixed. Equations [7] can be used as a check for the absence of errors in the computa- 
tions. 

In the language of electrical power engineering the quantities A, are m-times the 
“symmetrical coordinates” of the quantities (L,u*)°’’. Previous authors® have expressed 
the different roots of a polynomial equation in terms of their symmetrical coordinates, 
but these investigations have no bearing on the present problem. 

4. The computation of H(Z’) for symmetrically arranged poles Z; . Let the unit 
circle be divided into g equal sectors. g is either equal to m or a multiple of it. 


Let the radius of the unit circle be divided into 1/p) equal parts (1/po is an integer). 
In Fig. 1 this subdivision is shown for m = 8, g = 16, and 1/po = 4. The subdividing 
radii and circles determine “grid’’ points by their mutual intersection. Let equation 
[4] be rewritten as 


») - 3A: 
H(Z’) = Lz 7 (14) 
where 
Rag (v= 1,2,---,g) (15) 
and 
Al = af” + tal” = Ags , if v/k’ is integral, 
(16) 
Aj =0 if v/k’ is not integral. 





‘C. L. Fortescue, Method of symmetrical coordinates applied to the solution of polyphase networks, 
Trans. A.I.E.E. 37, 1315-1327 (1918). 

7C. F. Wagner and R. D. Evans, Symmetrical components, First Ed., McGraw-Hill, New York and 
London, Ch. XVI, 328-344 (1933). 

8C. L. Fortescue and G. Calabrese, L’applicazione delle coordinate simmetriche alla risoluzione delle 
equazioni algebriche, Atti del Congresso Internazionale dei Matematici. Bologna, p. 159 (1928). 
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Ao=A'g 
AsO © A;=0 
Az = Ae 6) 6) Ai = Ae 
A,=0 Os A =0 I 
Pp 
Py ec 
Ag? As O > ©) Ag* Aig me 
2 
A's 
] 
A'=0 
(e) oO ‘ 
A5=A 
me 7 ix 
Ag Aig 
Ay=0 — Ai3=0 


Let 


be a grid point, so that c is an 
ing 





Ag = Ai, 


Unit Crrc_e AND Grip PoInts 


Zz’ = pe” (17) 


integer. By substituting (15) and (17) in (14) and introduc- 


h=v- ce, if o> c| 
(18) 
h=gt+v-e, if v< cf 
we obtain after a few transformations 
H(Z’) =e ***"*H(Z"), (19) 
where 
ae) = } oe. (20) 
=i p —e 
Let the conjugate of H(Z’) be written as 
H*(Z’) = —(p + ir). (21) 


Then (—p) is the radial and 


field at the grid point Z’. The 


(—7r) the tangential component of the intensity of the 
use of [16] and a transformation of [20] results in 


h=9 h=g 
p= dafirr, — Dd. afs?t,, (22) 
h=1 h=l1 
h=9 h=@ 
oe? offln + 2 ion, , (23) 
A=! A=1 
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where 


and 
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th 


~ 1+ p” — 2p cos (2rh/g)’ 
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cos (2rh/g) — p__ (24) 
1+ p — 2p cos (2rh/g) 


sin (2nh/g) (25) 


Equations (22) and (23) are linear forms of the type of equation (8), because the values 
of r, and ¢, can be tabulated once and for all. Tables 1 and 2 are examples of such tables 


Values of r;, for 1/p 


TABLE 1 
8andg = 16. 








0'375 | p 




















































































































2h /g p =0 p = 0'125| p = 0'250| p = = 0/500 | p = 0'625| p = 0'750| p = 0’875 | p = 1'000 
0 +1 00 000) + 14 285 +133 333 +1'59 998 +2’00 000 +2’66 657) +4’'00 000) +8’00 000 * 
22°30 - 0’92 388) + vor 814 +112 210 41/22 597|+1'29 981 +126 768| +0’98 421) +0’32 843) —0'50 000 
; 45 +0 70 711 +0769 394 +0'64 477 40°54 417 +0'38 150) +0'16 204|—0’08 547) —0’31 736) —0’50 000 
67°30 +0’'38 268 +0'28 010 40'15 230 +000 900| —0’13 527|—0’26 562 0/37 160] —0'44 923, —0’50 000 
90 0’00 000 -0'12 308 023 529 032 877 —0’40 000 —0'44 944|—0'48 000] —0’49 357 —0’'50 000 
12°30’ || -0'38 268] —0'45 684|—0'50 460|—0'53 073|—054 064|—0'53 916| -0'53 015|—0'51 644| -0'50 000 
135 _0'70 71 060 784 0°67 590| —0’64 760) —0’61 680) —0’58 567) —0’55 547| —0'52 683) —0’50 000 
157°30 ; —0’'92 388 _0'84 140 -0'76 904 0°70 840|—0’65 500|—0’60 848|—0’56 774|—0'53 183) —0’50 000 
180 —1'00 000 _0'88 880 ~0'80 ‘000 _0'72 727) —0’'66 667) —0’'61 538) --0’57 143, —0’'53 333) —0’50 000 
—— Se. —— | — | 
*The value of this limit depends on the direction of approach. 
TABLE 2 
Values of ty for 1/pop = 8 andg = 16. 
= secs: = : —— 
Qn} = 0 p = 0'125| p = 0'250| p = 0'375| p = 0'500| p = 0'625 | p = 0'750 p = 0'875| p = 1/000 
0 ; 0 = 0 - a 0 —_ 0 0 | * - 
22 30” 0’38 268 0’48 771. 0’63 722 | 0’85 475 | 1117 347 | 162 311 2'16 607 2’57 126 | 2’51 383 
45° 0’70 711 0’84 295 | 0°99 741 | 1/15 863 | 130 252 139 541 140 903 | 133 874 1’20 713 
67°30’ 002 388 00 427° 1'06 053° “Vos 232 | 106 522 | 101 273 093 466 0’84 301 0’74 831 
. 7 ES eee | | 
90 1/00 000 | 098 461 | 094 118 | 0’87 671 0’80 000 | 0’71 910 0’64 000 0'56 637 0’50 000 
; . es a ial = 
112°30’ 0’92 - > = Anticadlhavendl 0'64 714 | 056 587 | 049 433 0’43 242 | 0’37 937 | 0’33 409 
135 0’70 711 059 301 | 0’49 935 | 0’42 318 0’36 131 0’31 O88 0’26 954 | 0’23 546 | 0’20 711 
157 30” ; 038 “268 0730 698 | 0/25 070 | 0’20 &71 | 0'17 608 | 0’'15 034 0’12 980 | 0’11 314 | 009 946 
wr || o | o | o | o | o | 0 ; o | o | o 


value of this limit 


| 
| 





depends on the direction of approach. 
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for 1/p.) = 8 and g = 16. Since r, = r,_, and t, = —4t,_, , it suffices to tabulate the 
values between h = 0 and h = g/2. 

If greater accuracy in the location of the roots is desired than determined by the 
mesh of radii and concentric circles, interpolation methods may be employed. H(Z’) 
is analytical and (dH(Z'))/dZ’ ~ 0 at a root point, unless (dG(uZ’))/dZ’ = 0 which 
corresponds to a multiple root. Consequently, in general, complex linear interpolation 
for H*(Z’) is admissible, if the values of H*(Z’) at two grid points are known. Other 
possibilities are repeating the described procedure for different values of u, or for equa- 
tions whose roots are shifted with respect to the roots of the original equation, or for 
equations whose roots are powers of the roots of the original equation, etc. 

After this paper had been completed a numerical method became known’ that uses 
Lucas’ principle to find with great accuracy “the roots of polynomial equations with 
exact coefficients when one has a first approximation as a starting point.’’ Since Salzer’s 
approach is entirely different from that in this paper, a combination of the two methods 
in question is not possible without a special investigation into its practicability. 

5. Polynomial equations with real coefficients. It will be shown in this Section that 
real values of the residues A, for a symmetrical arrangement of the poles along the 
unit circle can be obtained with the aid of a well-known conformal transformation, if 
the coefficients of the polynomial (1) are real. For real values A, equations (22) and 
(23) simplify considerably. 

By the relation 


(W + 1)(Z’ +12) = -—2 (26) 
the upper half of the Z’-plane is conformally transformed into the interior of the unit 
circle of the W-plane. The point Z’ = © corresponds to W.. = —i. The field configura- 


tion given by equation [4] in the Z’-plane, correlates to the field configuration given 


by the equation 


" F 7 s=m __ 7 — s=m " 
H'(W) = DX w-W.~ We > A, (27) 


in the W-plane, if the poles W, in the W-plane correspond to the poles Z{ in the Z- 
plane according to [26]. Saddle points of the potential correlate in the two field maps, 
hence we obtain the condition 
H'(W,) = 0 (28) 
if W, and Z! correlate. Consequently, the roots Z} of the equation G(uZ/) = 0 can be 
located by a field exploration according to Section IV carried out in the W-plane. 
The steps described in this Section so far can be comprised in a single procedure. 
Let the uniform distribution of the points W, including W. = —7’ be given by 


W, = an fg! (Be 0sa+8) (s _ a 2. ooo a 2). (29) 


) 


From (26) we obtain 


, — ___ Sin [2ns/(n + 2)] _ a ne 
4. = 1 — cos [2ms/(n + 2)]’ = 1,3, ate (30) 


%Herbert E. Salzer, On calculating the zeros of polynomials by the method of Lucas, J. Research Nat. 
Bur. Stands., 49, 133-4 (1952). 
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The substitution of (1) in (5) gives 


A,= bs (u*L,)-*h, ’ (s = 1, 2, “oe oe > 1), (31) 
where 
Ze (s = 1,2, --- ,n +1) 
h, =k BZ)’ (32) 


(q = 0, 1, 2, -++ yn) 


A, is a linear form of the type of equation [8], because the quantities *h, can be com- 









































TABLE 3 
Values of hg forn = 6 (k = 8). 
s=1 oa 2 | s=3 Pome s=5 eee | s=7 
“q@=0 || 4002 siz6| 1 | +4'97 487 -8 +4'97 487 -1 | 4002 512 6 
“q@=1|| -006 060] +41 | -200 006 | +2'06 066 -1 | +006 066 0 
q=2 || +0714 645 ~1 | +4085 355 | 0 40/85 355 at i +0'14 645 
q=3 || -o35 355 | +H | —0'35 355 24 0 +0'35 355 -1 | +035 355 
“q=4 || 4085 355 || =| 40/14 645 | 0 +0'14 645 -1 | +4085 355 
o=5 || -206 066 | +1 | 0706 066 0 | 0 +0'06 066 0 -1 | +2706 066 
q@=6 || +407 488 | ma 4 +0'02 512 e | 0 | +002 512 6 -1 | +497 488 





puted and tabulated once and for all. Table 3 contains the values ‘h, for n = 6. The 


residue at W.. = W,42 = —7 is given by 
s=n+1 

Aus = — >, A, = —k""'L,. (33) 
s=1 


The double statement in (33) can be used as a check for the absence of errors in the 
computations. 

6. Methods to obtain real residues A, for equations with complex coefficients. It 
may be sometimes desirable to work with real residues and to pay for this advantage 
with a larger number of poles. Two methods are suggested which achieve this and in 
which, at the same time, all roots or numbers closely associated with them are within 
the unit circle. 

Method 1: The equation 

t=2n q=n hen 
> kz = | = Lz|| 3 Laz | = 0 (34) 
t=0 e=0 = h=0 


has real coefficients and its roots are the roots of the original equation and their con- 


jugates. Equation (34) can be dealt with according to Section 5. 


Method 2: If G(Z) happens to be a polynomial of even degree in which any pair 
of coefficients equidistant from the two ends are complex conjugates, the residues A, 











326 L. TASN Y-TSCHIASSN Y 












[Vol. XI, No. 3 


computed from (12) turn out to be real. If the star denotes the conjugate, a subsidiary 
polynomial complying with this condition is 


> N,2" = | = wig2" |. ¥ (u"*L* a" (35) 


q@=0 h=0 
The (2n) zeros of the subsidiary polynomial are the n zeros of the original polynomial! 
G(uZ’) and its n conjugate reciprocals. 
If the values N, are computed and substituted in (12) for m = 2n + 2, we obtain 
after some transformations 


: .— Qrsq |” 
(—])’A, = | & on L.) COs 5 aa | 


[ h=n ) 
a 2 

h . “a7sq 4 ' eal 

T | 2 [,) sin > | , (s = 1,2, --- , 2n + 2). (36) 


i+ 2 
Equation (36) which is very closely related to linear forms of the type of equation (8) 
permits the direct computation of the values A, from the coefficients of the original 


equation. A check for the absence of errors is the first statement of equations (7). 











A METHOD OF SOLUTION OF THE EQUATIONS OF CLASSICAL 
GAS-DYNAMICS USING EINSTEIN’S EQUATIONS* 


BY 
G. C. McVITTIE 


University of Illinois Observatory, Urbana, Ill. 


Summary. It is known that Einstein’s equations in general relativity provide 
explicit expressions for the density, pressure and velocity of a perfect gas in terms of 
the coefficients of the metric (the potentials) and hence in terms of the coordinates. 
Using orthogonal space-times, the expressions involve four potentials only between 
which consistency relations hold. It is shown how degeneration of the Einstein equations 
to Newtonian hydrodynamics provides general solutions of the equations of classical 
gas-dynamics for motions which may be either of constant or of variable entropy. 
The consistency relations are obtained in the general case. As an illustration, one- 
dimensional gas-dynamics are discussed and it is shown how the consistency relations 
are manipulated. The solution in which one or other of the Riemann variables is constant 
is obtained as a special case and motions of variable entropy are also attained. 

1. Introduction. It is well-known that Einstein’s theory of general relativity is a 
generalization of the Newtonian mechanics of a continuous fluid but, as far as the 
present author is aware, it has not hitherto been realised that Einstein’s theory can 
serve as a tool in classical gas-dynamics. The object of the present paper is to show 
how this comes about. The solutions of the equations governing the motion of a gas 
which are obtained impose no limitations on the magnitude of the gas-velocity nor do 
they pre-suppose that adiabatic conditions prevail. But they do imply that the gas is 
perfect and non-viscous and that its motion takes place under the influence of its pressure- 
gradient alone. 

The Newtonian absolute time will be denoted by 7 and the rectangular coordinates 
in Newtonian absolute space by (X, , X. , X;). The three velocity components of the 
gas will be written (U, , U. , Us). The summation convention will be used throughout, 
repeated greek indices running through the values 1 to 4, whilst repeated latin indices 
will take the values 1 to 3. The letters J, m, n will stand for any cyclic permutation of 
the numbers 1, 2, 3. With these conventions, the classical equations of motion of a gas 
are 

al 
o7 


ov, Ss. ee eee ee 


te Me 3 


where (F, , F, , F;) is the force per unit mass of gas, excluding the pressure-gradient 
force; and the equation of continuity is 
Op 
oT 


The equations of motion can, with the aid of the continuity equation, also be written as 


Op U; _ 
+ Fy = 0. (1.01) 


0 


: 0 5 ae sr 
ap (el i) + ax, (pU,U; + 6p) = pFi, (¢ = 1, 2, 3) (1.02) 
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where 6,; is the Kroenecker delta. The four equations (1.01) and (1.02) express the 
conservation of momentum and of mass in classical gas-dynamics. The thermodynamical 
quantities usually associated with a gas-motion are the velocity of sound, a, and the 
entropy of unit mass of gas, which may be defined as follows: If the specific heats at 
constant pressure and at constant volume are c, and c, , respectively, and their ratio, 
c,/c, , is denoted by k, then 


a=k® 
p 


(1.03) 


’ 


The entropy, S, of unit mass of a perfect gas, whose temperature is @ and whose equation 
of state is p = Rpé, is most conveniently defined thus: Let dQ/dT be the rate of change 
“following the motion”’ of the heat-content of unit mass of gas, then 


IB = 3,5 + po (*), (1.04) 
whence, using the equation of state and the relations 

R= Jc, —c,) = (k — le,, 
it follows that 


fa a eee a 
paT ap (108 P k log p). 


In thermodynamics, the differential of the entropy is defined by dQ/6, but, for our 
purposes, it is more convenient to write 


dS = dQ/(c,@). 


Hence 
dS d r 7 - 
aT = aT {log p — k log p}, (1.05) 
or 
S = log «x + log p — k log p, (1.06) 


where « is a constant, characteristic of each separate unit mass of gas. In particular, 
if all unit masses have the same constant value of «x, the pressure and density of the gas 
are related by the isentropic relation 


p = kp’. (1.07) 


It is now necessary to summarize briefly the notions of general relativity which will 
be required*, but the reader who is not interested in this aspect of the matter may 
pass directly to Equations (2.06) to (2.09) and content himself with the process of 
verification described in the text immediately following these equations. In general 
relativity, an event is specified by four coordinates (x*, x’, x’, x°) of which the first 
denotes the time, and the other three, the place, of the occurrence in question. The 


*For a more detailed treatment see, e.g. R. C. Tolman Relativity, Thermodynamics and Cosmology, 
Clarendon Press, Oxford, 1934. 
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events which, for example, constitute the history of the motion of a perfect gas are 
regarded as mapped in a four-dimensional Riemannian space-time whose metric is 
ds’ = g,, dx” dx’. 
The velocity of the gas is represented by a four-dimensional vector (u*, u’, u’, u*) satisfy- 
ing the equation 
1 = g,yu'u'’, (1.08) 
whilst the mechanical quantities that are the counterparts of those whose partial de- 


rivatives appear in the equations (1.01) and (1.02), are the components of the energy- 
tensor 


Tt” = pv — o’ 5, (1.09) 


f 


+) 


where p, p are the (invariant) density and pressure, g*” are the contravariant components 
of the metrical tensor g,, , and c is the velocity of light. The energy-tensor and the 
metrical tensor are connected by the ten Einstein equations, viz. 


—8r7T" = G — 38°C (1.10) 


where y is the constant of gravitation, T? = g,,7,,, G? is the contracted Riemann-Chris- 
toffel tensor and G, the invariant curvature. The last two tensors can be expressed in 
terms of the g,, and their first and second derivatives with respect to the coordinates, by 
the rules of the tensor calculus. 

2. Newtonian approximation to Einstein’s equations. It is sufficient for our purpose 
to consider orthogonal space-times whose coefficients differ only slightly from those of 
the space-time of special relativity where 


1 
Gu = 1, Gis = Ges * Osa = ~ Juv = 0 (u # v). 


We write « = 2y/c’ and assume throughout the calculations that terms of order ¢ higher 
than the first are negligible. An orthogonal space-time of the required kind has metric 


ds’ = D(dzx‘)? — A(dz')? — B(dx’)? — C(dz*)’ (2.01) 


; {1 + Arde + ‘\h 
ras les} 
cn fit mle) 


and 9, ¢; , ¢2 , ¢3 are functions of all four coordinates (x* , x’, x’, x*). The right-hand 
sides of (1.10) have been calculated by Dingle* for general values of D, A, B, C. If terms 


where 
D=1 — 4reg, 


T 
iy 


2 


i) 


(2.02) 


a, | 





*H. Dingle, Nat. Acad. Sci., Proc. 19, 559-563 (1933). Also given in R. C. Tolman, loc. cit., pp. 
253-257. 
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of order ¢ are alone retained in his formulae, the following approximate forms of Ein- 


stein’s equations are obtained: 


lI 
- 
| 


4reT™ 


7 9” 
treT* = 4re se be + : en T od}, 





(2.03) 
a2 42 2 \ 
S11 oO” Dn Dn O Om OD» 
AreT” = 4e\ —7 2 (6 + r an ) 2 < 2 x 2(? 
(Ox) c (Ox ) (dx”)* 
bre Jo: é l > ( a a 3 | 
Te — —476€) 2) 3 — “a 
| ¥ C ad (ax”)? (d2") Pils 
where 
‘ a” 0° o° 
V* = 3 ——ang ——z, 
(dx')” + (Ox")” + (dx°)” 
Using also (1.09), the preceding equations become 
e OG, ] 
pu u —_—_ —_—- , vs 
Ox O2 
; 
p . = om ( Ym a es) 
1 = ¢ a < 
dx* da' | 
f (2.04) 
a a” Om T+ On O° Om 0°¢, 
pu) +p = -2 a (6 +* ts) +a + | 
(Ox) Cc (Ox ) (dx")" | 
| 
4\2 P 2 l f oO” a ) | 
(Ww ) — > = —Vo-—- — mee fe 
wre e . e a vax”)? +t (Ox ) jis J 
where, to a sufficient approximation, the velocity four-vector satisfies 
P 1 f ae a ul = 
(u’)? —- = V(u)" + (uw) + (uf = ]. (2.05) 
Cc \ 


There are eleven equations in the set (2.04) and (2.05), whose left-hand sides contain 
the six functions of the coordinates p, p, u‘, u’, u’, u’; whilst their right-hand sides 
involve the four functions ¢, ¢; , ¢2 , ¢3 . Obviously therefore there must be additional 
relationships between these two sets of functions and, in principle at least, it is possible to 
eliminate p, p, u*, u', u”, u® from (2.04) and (2.05) and thus obtain differential equations 
involving only ¢, ¢: , ¢2 , and gy; . Such equations will be called consistency relations 
because, unless 9, ¢; , ¢2 , 3 Satisfy them, the eleven equations (2.04) and (2.05) will 
be mutually inconsistent. We shall show later how these consistency relations are ob- 
tained and how they are manipulated in the calculations: we must first derive from 
(2.04) and (2.05) the corresponding equations in Newtonian hydrodynamics. This is 
done, as usual in relativity theory, by neglecting terms of order 1/c’ and, to this end, 
it will be assumed that ¢, 9, , ¢2 , ¢3 , and the four u* contain no terms of order c’. The 
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Newtonian velocity components (U, , U. , U;) will be the degenerate forms of the 
ratios (dx'/dx*, dx?/dx*, dx*/dx*) = (u'/ /ut whe. u®/u*) when terms of order 1/c’ are 
neglected; and the coordinates (z*, 2’, ?. '2°) will become the Newtonian coordinates 
(T,X, , X_, X,). Neglecting terms of order 1/c’ equation (2.05) reduces to the statement 


that u* = 1 whilst the following set of equations is obtained from (2.04): 

, 2 09, a i 
pU,.U,, = ax, aX,’ (2.06) 
oo. en (2.07) 

ei ee , 

‘2 _ Fe , Hen , Hm 

pl l + Pp — oT” + ax? ox?’ (2.08) 
p= —V'¢. (2.09) 


In these ten equations p, p now stand for the Newtonian density and pressure respec- 
tively. It is easy to show by direct partial differentiation that, if the quantities given 
by (2.06) to (2.09) be substituted into (1.01) and (1.02), these equations are satisfied 
identically provided that F; = 0, (¢ = 1, 2,3). Thus a solution of the equations of classical 
gas-dynamics has been obtained for the case when the gas is moving under the action 
of its pressure-gradient alone. The solution however involves the four functions ¢, ¢,(i = 
1, 2, 3) of the coordinates X; and of 7, which are not independent but are subject to 
consistency relations that may be obtained as follows: 

Elimination of p and the three U; from equations (2.06), (2.07) and (2.09) yields 


the three equations 


Te = ( xm fe) / ve. 2 

aX, aX, \dX,, dT dX, aT bila CH 
Again elimination of p, p and the three U, from (2.07), (2.08) and (2.09) yields three 
equations, only two of which are however independent, 


( 3” 3 Je 4 oe (0: — ga) = {(52%e) - (_#* ) | ve. (2.11) 
axe axi/** axt™ ~~ ™™ lk, a ax, ; 


The ten equations (2.06) to (2.09) may then be replaced by the five consistency relations 
(2.10) and (2.11) together with the further five equations 


ee a¢ / 2 _— ] 
U; = ax. aT V'¢; (¢ = 1, 2, 3), 
| 





p= —V'¢, 


; | 
7  < g 3 - 1 3 0’o, 1 3 (# He) Fi | 
site av | 3V(Le ) 3 Laxit3 > ax, aT ve h 





where, in the formula for p, it is presumed that ¢, , ¢g2 , ¢; have been chosen so as to 
satisfy the consistency relations. It is worth noticing that p and 42g are connected by 
Poisson’s equation and therefore 4x¢g is the gravitational potential of the distribution 
of gas. The corresponding gravitational force has however been neglected in determining 
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the motion of the gas since, as we have seen, equations (2.06) to (2.09) imply that no 
force other than the pressure-gradient is acting. 

3. One-dimensional gas-dynamics. As an illustration of the foregoing general theory, 
we consider certain types of motion in which the velocity of the fluid is parallel to a 
given straight line and in which the pressure and density vary spatially only with respect 
to distance measured parallel to this line. If the line is chosen to be the X,-axis and 
distance along it be denoted by X, motions of this kind are defined by assuming that 
all the variables in equations (2.10), (2.11) and (2.12) are functions of X and 7 alone. 
It then follows at once from (2.12) that U, = 0, U; = 0, and that only U, = U survives. 
The consistency relations (2.10) are identically satisfied, whilst the relations (2.11) 
are also satisfied by taking 


8, _ ( a ¥/z (3.01) 
aX’ = \aX aT ax” 
Hence formulae (2.12) reduce to 
ST ae, xen d"y Fi dy 
i aX dT/ ax” 
a9 ; 
= =, (3.02 
p ax 3 J ) 
7 _ dy ( ay \/22 
P= —ap? tT \ax ar) / ax® 





which provide the general solution of the problem of one-dimensional gas-dynamics 
which we are seeking by the method of Einstein’s equations. The solution is given to 
within an arbitrary function ¢ of X and 7’, which, from the purely mathematical stand- 
point, can be chosen in any way we please. But it is important to notice that not every 
such selection will give a physically acceptable solution as may be seen by considering 
the choice 


g = A cos(X — qT), 


where A, g are constants, which leads through (3.02) to 
U=4q, p = A cos (X — qT), p = 0. 


Thus the velocity of the gas is constant, its pressure is zero and its density is alternately 
positive and negative, vanishing at places where X = q7’ + (n + })a, (n, an integer). 
Such a distribution of gas would clearly not be considered physically admissible. 

The method of solution of equations (3.02) therefore consists of finding, by trial 
and error, a function ¢ that will correspond to a physically significant situation. The 
process of selection may be guided by some a priori requirement respecting the mathe- 
matical forms of U, p or p. As an example, let it be required to find all permissible den- 
sities and pressures for a perfect gas moving in such a way that its velocity is given by 


U = (n+ 1)[X/T + ng/(n + 1], (3.03) 
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where is a pure number and q is a constant with the dimensions of velocity. Writing 
u = d¢/AX and substituting from (3.03) into the first of equations (3.02), it follows 
that u satisfies the first-order partial differential equation 


oe nq On 
% +n + n(% ate ;) % =0, 


0 . 
n= eo NG + ot} 


where / is an arbitrary function of the argument (X/T + q)T”. It then follows that 


g= ron( + ar + H(T), (3.04) 


whence 


where H is an arbitrary function of T and 


F = = [ x(F+ a) be me aX, 


7 being treated as a constant in this integration. Introducing a new variable ¢ by 
X - . 
f=\p + qji'', (3.05) 


and denoting differentiation with respect to 7 by a prime, equations (3.02) yield 


U = (n+ 1){¢T" — q/n+ DD}, 
= —7 (n+1) d’F 
ani ar”? (3.06) 


: ce . 1 F H”’'T7"-* 
Dp = n(n + 1)T {t a(F ") —_ Hr 


By means of these expressions for p and p, it is possible to calculate the velocity of 
sound in the gas and the rate of change of entropy of unit mass following its motion. 
Formula (1.03) gives 


™ payed (F\_ WT) /(_@P ‘ 
= e+ Dn" 7 rat n(n + }/( a), G8) 


whilst (1.05) with d/dT = 0/dT + U(d/dX) yields 





dS _ k(n + 1) 


dT T 


n-1Ji. d (E i= rr} /{ 2d (£) be Hr} 
+ T \s dg \% (n + 1)n(n — 1) . d¢\¢ (n + 1)n)° (3.08) 
Obviously, a physically acceptable case must have p > 0, i.e. F must be such that 
d°F/d¢* < 0 for all values of ¢ corresponding to the region of space occupied by the 


moving gas, and for all relevant times. Moreover the pressure p and the square of the 
velocity of sound must both be positive; hence F, H and n must be chosen so as to fulfill 
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these requirements. As an aid to the selection of F, H and n it is instructive to specialize 
formula (3.08) in two independent ways. Firstly, let it be assumed that H satisfies 


H”’ = (constant)7”™’, (3.09) 
where the constant may have the value zero. Equation (3.08) then reduces to 
lS k 1 l ( owl ] 
dS _ ku +1) +m - 1) edi 
dT 7 
Alternatively, let it be assumed that each unit mass of gas conserves its entropy as it 


moves, a condition expressed by dS/dT = 0 or 


{ J ihe = {ie d (Z)- 7 : {hf ” pe — ‘ 
{kin +1) + (n 1)}¢ de \e la & ie {k(n + 1)H” + TH’’’} 0. (3.11) 


If both (3.09) and (3.11) hold simultaneously, and F is arbitrary, n is determined in 
terms of k by 
n= —(k — 1)/(k+ 1), (3.12) 
and since for a real gas 
> 8 > i, (3.13) 
n must be a negative number. 
(i) Velocity of Sound a Linear Function of X/T. 


The effect of choosing particular mathematical expressions for F must now be in- 
vestigated and, as a first example, suppose that 


F=-A? (3.14) 
where A(> 0) and Xd are constants; and suppose also that (3.09) is satisfied by taking 
H = 0. Then (3.06) and (3.07) become 
U = (n+ 1I{tT” — g/n t+ )D} = t+ Y{X/T + n¢g/(n + DV}, 


p= AMA— DE MT? = AMA — IY(X/T + Qn, 

p= —A(A — In(n + NPT" = —AQA — Inn + 1(X/T + TO’, (3.15) 
- {ty + neg) {_k iG + I vy T+q 

a= )- mn § | =)- \™ yp (X, 1), 


whilst dS/dT is given by (3.10). Thus the choice (3.14) for F leads to an expression 
for the velocity of sound which, like that for U, is linear in X/T7. But p, p and a’ must 
all be positive in a physically acceptable solution and therefore 


n(n+1)<0 and A> 1. (3.16) 


The first of these conditions is clearly satisfied in the constant entropy case in which 
n is given by (3.12). If further, the motion is isentropic and the expressions (3.15) for 
p and p are substituted into (1.07), it follows, by considering the indices of ¢ and T 
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in the resulting formula, that 




































k(\ — 2) =X and kn+1)=1-n. 


The second condition is, of course, the same as (3.12) whilst the first determines \ in 
terms of k, viz. 


A = 2k/(k — 1). (3.17) 


Because of (3.13), \ is greater than unity as required by (3.16). Using (3.12), (3.17) 
and (3.15), the solution for isentropic motion is 


p- 2 (¥_#=1,) 
 k+i1\7 2 , 


met DA(K, 
‘ieee (k aoe 1)? T +q ’ 
(3.18) 
2A 


me Ee 


II 
aa, 
~ 
S| 
rn) 
ee” 
ho 
! 





k—1/X 

a=t=1(X 4 Q), 
In the conventional treatment of isentropic motion in gas-dynamics based on the method 
of Riemann*, the solution (3.18) is obtained by putting one of the Riemann variables 
r,s equal to a constant. But the Riemann method gives, in the first instance, explicit 
values of U and of a, whereas the present one simultaneously determines p, p, U and a. 
This is due to the use of the gravitational potential 4x¢ that corresponds to the density p. 
A second point of difference is that the isentropic condition (1.07) is introduced a priori 
in Riemann’s method and it is therefore not easy to modify the method when variable 
entropy motions are in question. Such cases arise, for example, in the motions of inter- 
stellar gas clouds which are losing energy by radiation, a problem that has been attempted 
by Burgerst using the classical treatment. Variable entropy motions however present 
no greater difficulty than do adiabatic motions if the method we are here presenting 
be employed. For example, the solution (3.15) with the conditions (3.16) corresponds 
in general to variable entropy, the rate of change of entropy following the motion falling 
off inversely with the time by (3.10). In such motions the rate of loss of internal heat- 


energy in ergs per sec. per cm’* of the gas is 
dE oP prom J dQ eee Rpe_ dS oe Pp ds 
dT aT = ~k—-1dT) 3 =0ok—1dT 


Using (3.10), and (3.15) there comes 


_ = , r 
ae = oN i 1) f—-n(n+ )}{kn+1)+a- »i(2 + a) er, (3.19) 


*B. Riemann, Oeuvres Mathématiques, Paris, p. 177, 1898. A summary of the method is given in 


G. C. MeVittie, Mon. Not. Roy. Astron. Soc., London 110, 224-237, (1950). 
tJ. M. Burgers, K. Ned. Akad. v. Wet., 29, 600 (1946). 
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By hypothesis A > 0, and, by (3.13) and (3.16),2 >k >1,A > 1, and —1 < n < 0; 
hence all the constant factors in dE/dT are positive provided that k(n + 1) + (n — 1) 
is also positive. If the gas is monatomic, k = 5/3, and n must lie in the restricted range 
—1/4 < n < 0; if the gas is diatomic, k = 7/5, and n lies in —1/6 < n < 0; and so on. 
Thus by a suitable choice of n depending on the value of k, equations (3.15) give the 
motion of a gas which is losing heat-energy per unit volume at the rate given by (3.19). 


(it) Velocity of Sound a Quadratic Function of X/T. 
As a second illustrative example in the choice of F and H, we consider 


F = —Ag log (¢ + B), H"” = 4ABn(n + 1)T"", . = 0, 


where A(> 0) and B are constants. Equations (3.06) and (3.07) then become 


U = (n+ 1){X/T + ng/(n + 1}, 


p= a (X/T + q + 2BT7")/(X/T + q+ BT")’, 


p = —n(n + 1) a (X/T + ¢ + 2BT")?/(X/T + q+ BT"), 


a’ = —n(n + 1)kK(X/T + g + 2BT")(X/T + q + BT”), 


where n = —(k — 1)/(k + 1). Since n < 0 for any real gas, it is evident that, as T 
increases, U and a will ultimately both be linear functions of (X/T + q), as in (3.18), but 
the formulae for p and p will not become functions of (X/7T + gq) alone. 

4. Further developments and conclusions. The one-dimensional gas-dynamics dis- 
cussed in the preceding section do not exhaust the applications of the equations (2.10), 
(2.11) and (2.12). The case of spherical masses of gas in motion has also proved tract- 
able and, by proceeding to the second approximation to Einstein’s equations through 
the inclusion of terms in ¢’, it has been possible to take account of the gravitational 
self-attraction of the mass of gas. The initial coordinate system employed in (2.01) is 
moreover not unique and an analysis of the different permissible coordinate systems 
has thrown light on the meaning of coordinate systems in general relativity and on 
the problem of gravitational waves in that theory. These and other related problems 
will be discussed by the author in a forthcoming publication. In the meantime, it has 
seemed appropriate to give the essence of the method in the hope that it will be found 
useful by other investigators. 
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LARGE DEFLECTIONS OF A CANTILEVER BEAM WITH 
UNIFORMLY DISTRIBUTED LOAD* 


By F. VIRGINIA ROHDE (University of Florida) 





In this note only beams of uniform cross section are considered. The basic assump- 
tions are that the deformations are elastic and that the bending does not alter the length 
of the beam. 

Barten [1, 2] and Bisshopp and Drucker [3] have obtained results for a concentrated 
load. In the case of the uniformly distributed load, the Bernoulli-Euler equation, 


dé M 


ds EI’ 





1 
p 
gives, with the origin at the free end of the beam, 


dM 7a 6 
ds ws cos 6= EI as 





This equation does not lend itself to any simple solution. Put 





. dé — ex d’6 = a 
6= Das"; == )na,s"'; 3 = > n(n — 1)a,3"”. 
nen® ds a=} ds a=? 
The boundary conditions s = 0, 6 = a, d@/ds = 0 give a, = a; a, = 0. Then 
dé ws ’ ; ‘ - . 
3 = — py (608 & cos 7 — sinasin 7), T= > a,s". 
as J) n=2 


Expanding cos JT and sin T and comparing coefficients gives 


Goss *@ Qarvs = GD, n= 0, 1,2, ---; a, = 0; 


a, = —w cosa/6EI; a, = a;,w sin a/30ET; 


etc. Thus 
@=art > as. 
Since dy/ds = sin 8, 
y= [sin 6ds = T, sina + T; cosa; (1) 


T, = 8s — ajs’/14 — azags'’/10 — --- ; 
’ , 7. 3 sa 0s 
T, = a3s*/4 + ags’/7 + (ay — a3/6)s'°/10 + -- 
Similarly, 


x = T, cosa — T, sina. (2) 


*Received September 23, 1952. 
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Finally, 
M = EI >> 3ka,,8""'; 
k=1 
a=—)a;L™. (3) 
k=1 
With s = L, a = wL’*/6EI, the first approximation to the maximum deflection is 


P ] , - 
y = ssina+ 4 a3s' cosa = wL*/8EI, 


which is in agreement with small deflection theory. 

For computation, various values of @ are chosen arbitrarily and corresponding 
values of wL*/EI are found by eqn. (3). Against these values are plotted the ratios 
x,/L, the horizontal projection of the deflection curve to the total length; y,/L, the 
maximum deflection to the total length; and y,/y, , the maximum deflection to the 
maximum deflection as found by elementary theory. The results are shown in the figure. 
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For wl*/EI < 2, values obtained by small deflection theory agree quite well with 
those obtained by large deflection theory. It is to be noted that the use of small deflec- 
tion theory amounts to increasing the factor of safety. 
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NOTE ON RECTANGULAR PLATES: 
DEFLECTION UNDER PYRAMIDAL LOAD* 


By WILHELM ORNSTEIN (Newark College of Engineering) 


Using the Euler-Fourier method, a direct procedure leading to the evaluation of 
Fourier coefficients in the computation of the deflection of thin plates is developed. 
Consider a thin rectangular plate of uniform thickness, placed horizontally on four 
supports and acted upon by an arbitrary distributed load (Fig. 1). To solve the well 
known differential equation of the rectangular plate, 
= V°V*w = plz, y), | (1) 





y —4 























ty 


Fia. 1. 


the deflection is represented in a form of a double infinite series whose every term satisfies 


the boundary conditions: 


. Zz . J . 
w= > > A, sin me 7 Sin ne 4 (2) 
Substituting this expression in the plate equation (1), we obtain 
2 2\2 2 
m n . y l—yp 
Ane( 5 s.) sin mx — sin nx - = —az,7— P(A, Y)- 3) 
ed Ang +B nner, ser MY ( 


Multiplying both sides of equation (3) by sin (m’rx/a)-dzx and integrating from 0 to a 
and then multiplying both sides of the equation by sin (n’xy/b)-dy and integrating 


*Received October 23, 1952. 
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from 0 to b we obtain 


/ 2 2 2\-2 ab na 

Aun = ae (™, + n.) i il p(x, y) sin mr : az | sin nr dy. (4) 
The expression (4) has been developed for any arbitrary, continuously distributed load, 
but it may, with small modifications, be used for a concentrated load. 

Application to pyramidal distribution of load. The lcad is symmetrical with respect 
to two central axes AC and BD, see Fig. 1. If loads of equal magnitude are applied 
at two points equidistant from the axis BD, at P,(z, , y,) and at P,(x., y,) then it is 
obvious that due to symmetry zx, + x. = a and 

sin (mrx,/a) = sin mr(1 — 2,/a) = sin (mrz,/a) (5) 
and similarly sin (nry,/b) = sin (nry./b). Starting with the determination of A,,, 
due to the loads OFH and OFG, we see that when the contribution of the partial load 
OAH is known, then it is only necessary to multiply that value by four to have the 
expression for A,,,, due to loads OFH and OFG. By the same reasoning, when we multiply 
by four the contribution made by load OHB to A,,, , the value of A,,, due to loads 
OHG and OFF is obtained. The maximum load on the plate is at the center, where its 
value is P. At any point (x,y) on OAH the load is 


9P» 
oe Se (6) 


i.e., a function of x alone. Substituting this value of p in the expression (4), we obtain 


2 


8P(1 — pw’) [m? , n*\"? re? | f°” 4 i x - 
Ann. = ‘ (= + * | | | sin nr dy | —sin mr : dz. (7) 
ny } 


nr Elab a b Jo oe a 


In the second integral of the expression (7) the lower limit yo is equal to yo 
which is the equation of the diagonal HF. Because of symmetry of the load, both m 
and n are odd numbers, so that after integration the expression (7) takes the form 


8P(1 — ww) (m , n’\? f'?z . x z Jz 
= —— u) —=+5 | = sin ma — cos nx — d\—}. (8) 
x EIn a b Jo @ a asa 


Setting rz/a = u, the integral of this expression becomes 





= bz/a 





Mee 


1 x/2 
> [ u sin (mu) cos (nu) du 
WT Jo 


x/2 


x/2 . 2 
= . [ usin (m + nu du + ae | usin (m — n)u du. (9) 
2r 2r 


/0 


“a /0 


Integrating by parts, we obtain the value of A,,,, for the load OAH 


+ « =e (# ni) | co Mn + be cos [(m — *)e/2 (10) 
= EI c” = n(m + n) n(m — n) ’ 


Interchanging a with b and m with n, the value of A,,, for the load OHB is obtained: 


4. = —2PU = #’) (% ny |= [(m + n)x/2] , cos [(n — e/2 (11) 
— a ' m(m + n) m(n — m) 





nr EI 
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Adding to the expression (10) the contributions made by the loads OAE, OCG and 
OCF, and further, adding to the expression (11) the influence of loads OBG, ODE and 
ODF, the term A,,,, for the whole plate is 


a. ee) Sass ae 2 2\-2 
Ring 0 ll (™ + ".) cos (mar /2) cos (an/0). (12) 
a EImn a b 


The expression (11) gives a trivial value of A,,, , because with both m and n odd, 
An» Would vanish, which of course is impossible. It is seen that the equations (10) and 
(11) are true only for m # n, which, however, does not give any practical result. More- 
over, we cannot set m = n, in the expression (10) and (11) because a value A,,,, equal to 
infinity would result. It is clear that the integration performed is true only when m 
equals n; hence we must go back to the expression (9) and set there m = n. With this 
substitution expression (9) will yield the following value 


1 w/2 ] 
a2 [ usin 2mu du = ys (13) 


For the load OAH, the term A,,, from (8) becomes 


1  .? -2 
An. = 4A=#) (1, + 1) (14) 
xm EI \a b 
and for the whole plate 
_ 10P0=#) (1, 4 1)” 
Ris — x mE! a’ + b? ° (15) 


Finally the expression for the deflection is 


16P(1 — yp’) (? | y{ . 2... 9 1 . Sex . Say } 
ee et FTP res eS ee ee eo 6 
u — r + 52) \sin | sin > + 3° sin — sin > + (16) 


Thus, through this operation a double series for the deflection of the rectangular plate 
under pyramidal load is reduced to a result involving but a single series. 


A RANDOM WALK RELATED TO THE CAPACITANCE 
OF THE CIRCULAR PLATE CONDENSER* 


By E. REICH (The RAND Corporation, Santa Monica, Cal.) 


Abstract. It is shown that the solution of Love’s equation for the capacitance of the 
circular plate condenser can be expressed in terms of the mean duration of a certain 
one-dimensional random walk with absorbing barriers. The interpretation as a random 
walk makes it possible to confirm the fact that the actual capacitance of the condenser 
is always larger than the value given by the standard approximation for small separa- 
tions, and yields an upper bound as well. In addition to its theoretical interest, the 


*Received November 21, 1952. 
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random walk appears to provide a practical means for the calculation of the capacitance 
by a Monte Carlo technique. 

1. Introduction. The purpose of this note is to point out an equivalence between the 
mean number of steps of a random walk, and the capacitance C of a condenser consisting 
of two equal, infinitely thin circular conducting disks of radius a, separated by a distance 
Ka. 
It is a classical result that for «x — 0, C ~ a/4x, and it is also known (Jeans [1]) 
that lim,...C = a/x. A more precise formula for small « dates back to Kirchhoff [3], 
who also discusses a paper by Clausius on this subject. A formula of Nicholson’s [5] 
said to give C for all «x by means of a definite integral has turned out to be incorrect. 
(Love [4] has pointed out fallacies in Nicholson’s reasoning. Moreover, it can be shown 
by direct consideration of Nicholson’s integral that it does not even give the correct 
asymptotic formula as x — 0). Present-day knowledge regarding the value of C for 
general values of x appears to be restricted to the following elegant result obtained by 
Love. 


LemMaA 1. If f(t) is the solution of the integral equation 


; i te f(t) 
f(x) — —— —, dt = 1 —I/k<2r< : 
f(z) of ice ( csi 2¢ 1/x®) (1) 
then 
Ca. (2) 
T 
where 
Z nl/« 
w=<] fd dt. (3) 
- —1/« 


It will be shown that both f(é) and yu have very simple interpretations in terms of 
a random walk in the interval | x | < 1/x. The results are expressed by the following 
theorem. 

THEOREM 1. Let 6;(i = 1, 2, ---) be independent random variables, such that 
0 < 0; < 2x, and Prob {6; < A} = 1/2r A, (O < A < 2m). Consider the random walk 
starting at x e [—1/x,1/x] whose 7-th step equals tan 6, . We will say that absorption 
occurs at step k if the k-th step is the first step landing outside [—1/x,1/«]. Then f(x) 
is the expected number of steps to absorption. 

Coro.LuaRy. Let the walk start at a random point in [—1/x,1/x] (i.e., a point 
drawn from a distribution with flat density). Then yu is the expected number of steps to 
absorption. 


2. Derivation. It will be seen that Theorem 1 is a corollary of the following result: 


THEOREM 2*. If 
(a) K(z,t) is continuous and > O0O(a<zx<ba<t<)b), 


(b) there exists a constant a > O such that 


| K(z,t) dtc l—a for all ze[a,b], 


“a 


*See Wasow [6] for a different proof of this result under a less restrictive hypothesis. 
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then there exists a random walk starting at 2 whose mean number of steps to absorption 


is f(x), where 


ja) -— | K@osOd=1 @s<z<yd). (5) 


The dispersion of the number of steps will be finite. 
First we derive the rather obvious 


LemMaA 2. If and 7 are random variables, taking on non-negative integral values, 


and 
P, = Prob {& > k} < Q = Prob {n > k} (k = 0, 1, 2, -->) 
then 
¢) <(n) and ¢’) < (n’). 
Proof: 


(n) = DO MQe — Qees) > DL MQe = Queer) + + DQwer = ps Q. . 


Letting N — , we see that the term to the right of the inequality approaches (7), 
showing that (n) = >—% Q, , and consequently (¢) < (n). The inequality for the second 
moments can be proved in the same way. 


Proof of Theorem 2. For each z e [a,b] construct a continuous function S(z,t) such 
that 


(a) S(z,t) > 0 (—-wx< t<o) 
(b) [ S(z,t) dt = 1 


(c) S(z,t) = K(z,2) for a<xt<b. 


For fixed x, the function S(z,t) should be thought of as a probability density. The random 
walk is defined inductively as follows: We start at 2, and stop when first landing outside 
[a,b]. Having taken n steps, and assuming the nth step lands us at y e [a, b], the next 
step is taken to the point A, where 


Prob {A <u} = [ S(y,t) dt. 


Let £ = &(x) = number of steps to absorption if the random walk starts at z, and 
g p Pp 
P, = P,(x) = Prob {& > k}. If before a given step occurs, absorption has not yet taken 


place, then the probability of absorption as a result of the step is not smaller than a. 
Thus P, = 1, and P, < Q, , where Q, = (1 — a)**. By Lemma 2, (&) = f(x) < 1/a, 
and (¢?) < (2 — a@)/a’. Dispersion [] = ((— — <é))”) = () — (&)’ < (2 — @)/a’. Thus 
f(x) exists, and has a finite dispersion. 
Now to show that f(x) satisfies (5). 
Define f(x) = 0 for x g [a,b]. Then for all random walks whose first step carries to 
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y, the mean number of steps to absorption is 
fz«|y =1+ fy 
Therefore 


a) ao 


| f(x | yS(a,y) dy = S(x,y) dy + | S(y) S(2,y) dy, 
which reduces to 
ab 
f(z) = 1+ | f(y) K(a2,y) dy 
as was to be proved. 

Theorem 2 follows from the observation that the tangent of a random number uni- 
formly distributed in (0,27) has the probability density (1/m)/(1 + 2”). The a of Theorem 
2 can be taken as 

sl/« 
dx 
l—1/r | = (2/mr) Arctan x. 


Jewel +2 


Thus 


f(z) < z mere ’ and Cs - for any x. 
2 Arctan x 2 Arctan x 

3. Further remarks. It is interesting to note that checks on the asymptotic values 
of C as x — o, and x — 0 are easily possible. 

For x — ~, we have of course f(x) — 1, as can be seen either from the random walk, 
or the integral equation (1). Therefore n — 1, and C — a/rz. 

For x — 0 the formula C ~ a/4« can be obtained by using a result of Kac and Pollard 
[2] which also yields the conclusion that C > a/4«. Kac and Pollard consider a con- 
tinuous random motion on the z-axis, where if the particle is at point z> at time ( , 
its position at time ¢ + f is given by a Cauchy distribution with semi-interquartile 
range ¢. We now make the observation that if one observes such a continuous random 
motion at the discrete time points ¢ = 0, 1, 2, --- , and records the observed positions 
at these instants, the result is the generation of a discrete walk with exactly the properties 
defined in the hypothesis of Theorem 1. 

Let f*(x) be the mean time to absorption of the continuous random motion, if it 
starts at z. As a simple consequence of the above observation we have f(x) > f*(z). 
Therefore 

ol/« 


up> «/2 | f*(ax) dx. 


l/« 


Using the result of Kae and Pollard for £{7'(a,b,t)} ({2], page 383) we obtain 


Therefore 


y _ aL ak . gla 6 
om > $8 ["" pry de = f. 0 
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In other words, C is always larger than the value predicted by the usual asymptotic 
formula for small separations. For x — 0 it is heuristically evident that f(r) ~ f*(2), 
so that C ~ a/(4x), as expected. f 

A simple, plausible, “physical” proof of the inequality C > a/4x can be given.* 
C equals the ratio of charge per plate to potential difference between plates. The capacitor 
can be considered as being cut out of a pair of infinite parallel plates. Before the cut- 
away portions are removed the ratio of charge to potential difference is a/4x. When 
the cut-away portions are removed each particle of charge moves away from the axis, 
and therefore the field intensity at each point of the axis is decreased below its previous 
value. Since the potential difference is the integral of this field intensity along the axis, 
the argument shows that the potential difference is reduced, and the capacitance there- 
fore increased from a/4x. 

4. Practicability as a Monte Carlo method. If « = ((¢ — u)’) is the standard devia- 
tion of £, and n walks are taken, yielding &, , & , --- &, as the observed values of ~, then 
a/n'’* will be the standard deviation of the empirical mean, and ¢€ = o/yun'”’ will be 
approximately equal to the ratio of standard deviation to the mean. 

We will restrict ourselves to small x, as it is conjectured that the n required for 
achieving a given ¢ steadily decreases as x increases. For small x, ¢ < 2'”’/a~ (2'*x), 


un > w/4x. Therefore, asymptotically for small x, 
e < (8/n)'”’. (7) 


A numerical experiment, using x = 0.1, and n = 1000 was made, using automatic 
computing machinery. The results were an empirical mean value of 10.101, and an e of 
approximately 0.034. For n = 1000 the right-hand side of (7) yields « < 0.09. 

5. Acknowledgment. The author wishes to thank Mr. P. Swerling for his helpful 
suggestions during the writing of this paper. 
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ON THE STABILITY OF LAMINAR BOUNDARY LAYER FLOW* 


By SIN-I CHENG (Princeton University) 


The stability of two dimensional small disturbances in laminar boundary layer 
flow has been extensively studied on the assumption that boundary layer flows are 
essentially parallel flows. The direct effect of the local pressure gradient on the calculation 
of the stability limits for incompressible boundary layer flow has been shown to be 
negligible under the approximation (aR)~* < 1, if the local velocity profile is used in 
the stability calculation. However, the assumption that the vertical velocity in the 
boundary layer flow plays negligible role has not received careful attention. 

There is a qualitative argument that under the Prandtl boundary layer approxi- 
mation, the variation of the mean flow properties in the z-direction within a few wave 
lengths of the disturbance is of the order of R~’, which is negligibly small compared to 
unity. Therefore, the contributions of such terms due to x-gradients of pressure and 
temperature in the stability calculation can be neglected as higher order small quantities. 
This kind of argument should be investigated more closely so far as the vertical velocity 
component is concerned even though the vertical velocity component is a small quantity 
of the order of R~'. The vertical velocity component produces a momentum transfer and 
an energy transfer across the boundary layer where both the disturbance quantities 
and the mean flow properties vary rapidly. The net effect of the transport processes 
may thus be much larger than the magnitude of the small agent that produces the 
transport processes. While the vertical velocity component of the flow and the gradients 
of the flow properties in the z-direction are small quantities of the same order the net 
effect of the former in the stability calculation will be shown to be much more important 
than that of the latter. 

From the two dimensional forms of the equations of mass continuity, momentum 
and energy with variable viscosity and thermal conductivity coefficient and constant 
specific heats, the linearized system of partial differential equations for the amplitude 
functions of the periodic disturbances of the type exp [a(x — ct)] is obtained with two 
independent spatial variables x and y. The fast varying, frequency dependent part of 
the disturbances as functions of x is represented by the factor exp (tax). The decay or 
growth of the amplitudes of the disturbances in the z-direction is slow and in the first 
approximation, the x-gradients of the amplitude functions can be considered as inde- 
pendent of x. The system of disturbance equations can hence be considered as a set of 
ordinary differential equations for f, ¢, 7, r and 6 with y as the only independent variable. 
Thus, we have the following equations: 


Continuity: 


“haa ; r p’ rv’ 1 oO 
— Fo SF (of + wi) (1) 
ap p p pa ap Ox 


*Received December 15, 1952. This study was supported jointly by the Office of Naval Research 
(U.S.N.), Mechanics Branch, Mathematical Sciences Division, and by the Air Research and Develop- 
ment Command, U.S.A.F. under Contract No. N6-Ori-270, Task Order No. 6 Project Number NR- 
061-049. The present paper is an abbreviation of Report No. 211, Aeronautical Engineering Department, 


Princeton University. 
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First Momentum: 


apli(w — c)f + w’d] + plfw, + wf. + vf’) + (ww, + vw')r = sR (iar + 1, 


] 4 2 . nie Qe 1 2 ” 
+ p (3 M+ 3 mT + 2af, — af) + (# += 3 ua al + tag’) + uf 


9 
+ , om (3 + 3 Af. + taf) + 2 (5 — 1)T,ad’ + aT’ (¢, + tad) + rs" | 


T . H (: + : rw, + 2 (7 - be’ fa, + tad) + (v. + we" | 


P I dy, f 2 ” l 2 / 
T RdaT (3 + 3 rw, Niitatlibial (} + 3 rt (2) 


Second Momentum: 


a’ p[t w — co] + aplv’d + ve’) + pv.f + apwe, + (wo, + vv’)r 


= © {Efi + taf’ + agree + Zia’o, — a] 


R 
2 p> 2 (us 

eft oheoan 

e oO My 3 \u 

1 dy In (f’ + + i 7 +: 2 Nad 4 3 — 1)7T’(f. +2 ni 
+. R aT A - ad, ia’) 37) +3 T Vs a. J 
1 1 du, fa. ' v . (4 2 \y , 2 pes ) 7 

Par (v, + w’)(6, + 108) + 3 tatpe + 3 (r 1)w,6 

de | (4 2 \y (: 2 7 nr’ 
‘v2 + 9 — ——s5 
T RaT \"** is+35*} + Ag + 3 7/eel? — ap (3) 

Energy: 


apli(w — c)0 + T’o] + plwd, + v0’ + T.f] + (wT, + oT’)r 
= —(y — 1)[p(f, + taf + ad’) + (w, + v’)x] 
; vy — 1)M*. 2 we ; a 
+ My = is + 3 u ) atu. fz + taf) + av’¢’] 
- ¢ — 1) wae + v'(f, + iaf)] + 2(w’ + v,)(f’ + ad. + ia} 
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State: 
vy r 0 
=—-+ a, 
p p 7 
where + = [(du./dT)/(du,/dT)], all other symbols are adopted from reference 1. 
For the convenience in reducing the equations the quantities Z; are defined in a 
manner slightly different from those in reference 1 and following the well known pro- 


cedure, these variables Z; are expanded into series of « = (ah)'” as follows: 
ZL, = f=a, + ex, + 
thy = ef'=z, + ex,’ + 
a 7 e'@=2, +e; + 
Zs - ¢’ a ‘i’ a ex,” + 
Z, = 6 ry + «x; + 
eZee = 60’ = 2%. + Xa. + 
1 1 W ) a ' 
€ Zr ™e _ +er, 
M 
- 1’ ( (1) ’ 
LZ. = Ve t. +e + --- (6) 
The mean flow velocity components w and v, the mean density and temperature are 
also expanded into Taylor series about the critical point y, where w(y.) = c. The ex- 


pansion of v is written as: 


av ou 


| 0 1 | 0 
n | — (pw) dy -* | — (pw) dy 
pJo OX R p 0& 


ae [vo + v,-€n + v2-(en)*/2 + ++] 
1 | O(pw) O(pw). .o 

Yo = ty, — A 92/2 + > 
p O& ~ 0& 


| | Sew 
; eo a = 
Pe oF p 


are of the order of unity for both compressible and incompressible flow. Substituting all 
these expansions into the disturbance equations and collecting terms of the same power 
of ¢, one obtains for the first approximation of the order of unity or e’: 


where 


a 
| 


dr 


+iz, =0 
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These equations are identical with those given in reference 1, and solutions can be 


obtained therefrom. 
For the second approximation of the order of ¢« one obtains: 


dx 
de 


= ’ 
vw; Pe (0) 
= T nx mae r3 
- P 
u v 1) 
(inary 1. 4 —_ t+ 
I YHic 





( a , , 7 2 
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ies ( dT Dent; aT (w clg + Tix ) 
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-~ iene + (r: ee Sat vel 
Vy) 7 Pp. 


= — ‘[2 _— (a (In ws ) 1 | iw’ xy” + Tinx;”) + iw! 


xy” 


to I. 


p . aT 


) T a a + T' x: } 
c Pe 


__ eed E = (aie My 
Y Pe dT 


oe —leo oo — i (0) —_ (se. Mi w) , 
Pies nM? x; Qo(y — 1)M*wix aT 2T x. 


(9) 


The homogeneous parts of equations (9) are the same as equations (8), and the inhomo- 
geneous parts are known functions in terms of the solutions of equations (8) and mean 
flow properties. The second approximation can therefore be obtained by quadrature. 
It is observed that the vertical velocity component enters into both the momentum 








1 
i131 


of the order of ¢’. 


‘: (0) 
Since 23; 
1 

and 25 


or x{)) and 233 


and the energy relations in the system of equations (9), while the z-gradients do not at 
this approximation. Terms involving z-gradients enter only at the next approximation 
Therefore, the effect of the vertical velocity component is more critical 
than the effect of the gradients in the z-direction in the stability calculation. 
= (dx\°/dn) = 0 for j = 3, 4, 5 and 6, it is only the two viscous solutions 
’ that are dependent on the vertical velocity component. 
x are all v-dependent but they enter into the characteristic value problem as higher 
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order small quantities even in the compressible case. The major effect of the vertical 
velocity component in the determination of the stability boundary is through x;'’ and 
T52'. 
Suppose we take all z;; from the e series and consistently take all six solutions 
to the order of ¢ i.e. x;;’ + ex;;’ in the boundary value problem, the v-dependence of 
23; and z;,’ indicates an inconsistency of the simplification of assuming the boundary 
layer flow as parallel flow. Fortunately, all the previous investigators are satisfied with 
the first approximation of the two viscous solutions x;;’ and x3,’ while they used x;3’ + 
ex;,. and x;;’ + ex ;,’or the two equivalent inviscid solutions or the inviscid solutions 
corrected for viscosity in the boundary value problems. Hence, the stability boundary 
as determined by any of these methods is independent of v and their results are consistent 
with the assumption that boundary layer flows are parallel flows. Therefore, within the 
order of approximation attempted by previous investigators, the stability of the laminar 
boundary layer is determined only by the local flow properties for both the compressible 
and the incompressible flow. 

The accuracy of the quantitative determination of the stability boundary as carried 
out in references 2, 3 and 4, however, can not be improved merely by taking more terms 
in the ¢ series without including the effect of the vertical velocity component. It is 
unfortunate that, in some practical cases, the parameter e« may be only 0.1 near the 
minimum critical Reynolds number and as such the second approximation of the order 
of e should better be considered for accurate determination of the stability limit and the 
initial amplification rate, in which cases, the effect of the vertical velocity component 
must be included. 

In addition, at high Mach numbers, the vertical velocity component in the boundary 
layer is of the order of M*/R, and may enter the stability problem of the laminar bound- 
ary layer even in the first approximation. The stability of the hypersonic laminar 
boundary layer, therefore, requires careful investigation. 
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STRESS-STRAIN RELATIONS, UNIQUENESS AND VARIATIONAL THEOREMS 
FOR ELASTIC-PLASTIC MATERIALS WITH A SINGULAR YIELD SURFACE* 


By W. T. KOITER (Technical University, Delft, Holland) 


1. Plastic stress-strain relations. The state of stress at any point of a continuous 
medium is described by the stress tensor ¢;; and may be represented by a point in nine- 
dimensional stress space. It will be assumed that no yielding occurs if the stress point 


*Received Dec. 22, 1952 
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lies within a convex domain in stress space, which will be called the elastic domain. 
The boundary of the elastic domain is the yield surface which is called regular if it is 


described by an equation 
f(o:;) = 0, (1) 


where f is a regular (continuously differentiable) function, the yield function, of its 
nine variables, symmetrical with respect to o;; and o;; . The sign of f is chosen such 
that f(o;;) is negative in the elastic domain. For a so-called perfectly plastic material 
that yields under constant stresses, Prager [1] has shown under some very general 
assumptions that the plastic strain rate tensor is given by 


Ne iat of 
é;;/ =X ae,’ (2) 
where 
1 =O for f<0 ] 
ae 
and also for f =0, j= inn ai; < 0; t (3) 
, ‘ >: 
> 0— = ieee 5; = 
A > ( for f =0, f oa Ci; 0. | 





For a material with strain-hardening the relations (2) and (3) are replaced by 


3’ = 0 for {<0 
ui a se OF. :' 
and also for f =90, f= ac, "1 = 0; (4) 
ayant sg for f=0 pf=aAe,>o, 
OO ;; 00;; 


where h is a scalar function of stress, plastic strain and strain history. The plastic stress- 
strain relations (2), (3) and (4) may be called associated with the yield function. If the 
plastic strain rate is depicted in nine-dimensional stress space, these relations are ex- 
pressed by the geometrical statement that the direction of the plastic strain rate is 
given by the normal to the yield surface. 

In this paper singular yield surfaces will be considered of a type described by a 
number of regular yield functions f,(¢;;) (symmetrical with respect to o,; and a;;) 
such that the elastic domain is given by 


f,(o:;) <0 aaa l, 2, oe (5) 


and that yielding occurs as soon as at least one of the functions f, is zero. All points of 
the yield surface, where only one function f, = 0, are regular and the corresponding 
plastic strain rates are the same as for a completely regular yield surface f, = 0. Diffi- 
culties arise only on the intersection of two or more surfaces f, = 0 because for these 
stress states the normal to the yield surface becomes indeterminate. However, this 
ambiguity is largely removed if the criterion for unloading or yielding at such a singular 
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point is applied for each yield function separately. Physically this assumption is entirely 
plausible because each yield function represents a separate yield criterion. Its mathe- 
matical expression is for a perfectly plastic material 


‘ ~ of 
= Dr, (6 
€;; > L 0o;; i) 
where 
1 =0 for f, <0, 
, 5. Of, 
and also for f, = 0, j,a= a<& ™ 
. 00; (7) 
; : - Of, 
x = 0 for f, = 0, \ -—— ¢,, = @ J 
F ; O¢;; 


For a material with strain-hardening the relations (6) and (7) are replaced by 

— : oe 

é;; = Zz. a 7. (8) 

p=1 00; ; , 
where h, is a positive function of stress, plastic strain and strain history, and c, is zero 
if either f, < 0 or f; < 0 and unity for f, = 0 and f; > 0. 
2. The variational and uniqueness theorems. Current proofs of the variational and 

uniqueness theorems for a material with a regular yield surface and the associated 
plastic stress-strain relations are all of them based on the inequality (see [2]) 


CEG) H CigEii(2) — Wisse, > O, (9) 


where 635.1) , €:;<1) aNd o;;;2) , €;;2) are two arbitrary pairs of stress rates and associated 
strain rates at the same point of stress space. The equality sign in (9) holds only if 
0330) = 03:2) . This inequality is applied directly to prove the two’variational theorems, 
and the uniqueness theorem is shown to hold by means of the inequality, obtained by 


interchanging the suffixes 1 and 2 in (9) and adding 


ii €ii(2) — Ces€isay & O. (10) 


O87 yEsia -+- Oij(2)€84(2) ~~ Fij(i)€ii(2) al 


The proof of (9) is obtained by considering the elastic and plastic strain rates sepa- 
rately. Hooke’s law leads to the inequality 


C; ves) Hh Cems = 26 i 5 (1) €i5 (2 > 0, (11) 


where ¢;; is the elastic strain rate tensor; and the equality sign applies only if ¢;;,,, = 
04:2) - The similar inequality for the plastic strain rates 


o; 1 es + Cii(2 €is(2 — 265 (1) Eis (2 > 0 (12) 

is shown to hold by means of the stress-strain relations (2), (3) or (4). 
It will now be proved that the stress-strain relations (6)-(8) also lead to inequality 
(12) and hence ensure the validity of the uniqueness and variational theorems. For a 
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perfectly plastic material the left-hand member of (12) is calculated as follows 


; oe — re ong 
Ossa)eisia + O475(2)€i5(2) 26 5 5 (1) €55(2) 
n 
= | crete 22 + denotes 22 = Mp metien LE 
ad p(1yFij(a) a + p(2)%ij(2) a = Ap(29 Fis) a 
p=) Ci; C3; Ci; 


Il 


>» Dou toay + Apia F viz) ae rpc» f ou]. 
p=1 


The various possibilities for each term in this sum are listed below with the consequences, 
ensuing from the stress-strain relations (6), (7): 





fe <O—Asay = Ae) = O, (a) 
f, = 90, fou < 9, fey < OA, = Apes = 0, (b) 
f,=0 fruy = 0, fore») KOA DO, = Ava = *O, (c) f (13) 
f,=0 finn << fo = OA) = O, A) > 0, (d) 
f-=0, fia =0, fea =OAwW > 0, An > 9, (e) J 


It follows that only if (13d) applies a non-zero term occurs and this term is then always 
positive. The inequality (12) must therefore hold true and the uniqueness and varia- 
tional theorems are proved. 

The proof of (9) for a strain-hardening material with stress-strain relations (8) is 
entirely similar and may be omitted here. 

3. Application to Tresca’s yield criterion for a perfectly plastic material. In some 
problems, e.g. the thick-walled tube under internal pressure [3], Tresca’s yield criterion 
and its associated stress-strain relations afford a much simpler approach than von 
Mises’ yield condition and its associated Prandtl-Reuss relations. It will now be shown 
that Tresca’s yield criterion belongs to the singular yield criteria considered here. 

The most convenient formulation of Tresca’s yield criterion for the present purpose 
is that the shear stress component in any plane, defined by its unit normal n;*’, and in 
any direction n;”? in this plane (which is of course perpendicular to n{°’ shall not exceed 
the yield stress in pure shear. This criterion may be written in a form symmetrical 
with respect to o;; and o;; 


5 a.s{nin\® +nint?} —k <0, (14) 
where k is the yield stress in shear. Clearly the number of yield functions is now infinite 
but this does not affect the foregoing argument and the uniqueness and variational 
theorems remain valid. The plastic strain rates associated with the yield functions 
(14) follow from (6) 


R 1 , 
ef = 2. 5 Nas» fms n,” + nj?nt?}, (15) 
(a,b) 
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where },,) is only non-zero for those directions n;*’ and n;' for which the equality 
sign in (14) applies and moreover 

a: {nen + nin?) = 0. (16) 
It is obvious that the plastic strain rate, associated with a yield function (14), is a pure 
shear in the directions n;"’ and n; 

If the principal stresses are unequal, 7, > o2 > a3, there is only one pair of directions 
n\® . n\” for which the equality sign in (14) can hold; this occurs if ¢, — o,; = 2 k and 
the plastic strain rate is a pure shear with its principal axes along the major and minor 
principal stress axes. If two principal stresses are equal, e.g. ¢, > o2 = o3 , there isa 
one-parametric family of pairs of directions n;* , n;"’ , for which the equality sign in 
(14) holds. However, if the stress rates are such that oj = o; < o,, there is only one 
pair of directions in this family for which (16) holds true and the plastic strain rate is 
again a pure shear. On the other hand, if both o, — 2k = o. = o3; and a; = o; = a3, 
the plastic strain rate is indeterminate; the stress-strain relations only require that 
one of the principal axes of the plastic strain rate tensor coincides with the major principal 
axis of the stress tensor, that the two other principal plastic strain rates are both negative 
and that the plastic volume strain is zero. 

4. The slip theory of Batdorf and Budiansky. Considerable interest has been aroused 
by the slip theory of plasticity for strain-hardening materials, advanced by Batdorf 
and Budiansky some years ago [4], partly because it seemed to represent a new approach, 
entirely different from flow and deformation theories. However, this theory may be 
regarded as a special case of the present theory for materials with a singular yield surface. 
The number of regular yield functions, of which the yield condition is composed, is 
again infinite; the yield functions are similar to (14) 

] 
9 


— 


(< (b) 3 b) ) —_ 
o,;{n; nn; +n n; | — kas) < 0, (17) 


where k,,,, now represents either the initial yield stress in pure shear or (if the latter 
is higher) the highest previous value of the shear stress component with respect to the 
two mutually orthogonal directions n‘”’ and n;”’ . The plastic stress-strain relations of 
Batdorf and Budiansky may be expressed in the form 


. 1 (a)__ (b) 1 (b) (¢ (b (< (b) . 
¢;' = | ecntt ab) 4 inven? + nj ny tne ny” + ny n, }o,, IQs) , (18) 


where dQ,,,, is the measure of the three-dimensional set of pairs of mutually orthogonal 
directions n‘” , n‘” within an infinitesimal region, H,,,) is the characteristic strain 
hardening function of slip theory which depends only on k,,,, , and c{.,) plays the role of 
c, in (8). It is easily seen that (18) represents an example of (8) for the case of an infinite 
number of yield functions f, . Consequently the uniqueness and variational theorems 
are equally applicable to a material that follows the slip theory of plasticity. 
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SERIES EXPANSIONS OF RAYS IN ISOTROPIC, NON-HOMOGENEOUS MEDIA* 
By SAUL GORN (Aberdeen Proving Ground) 


The interest in air-to-air radio propagation and its possible use for high-accuracy 
distance measuring make it useful to have available recursive formulae for the coefficients 
of the power series expansions of rays. The following note has been extracted from a 
government report (see Gorn, [6]) with this purpose in mind. 

It will be convenient to use a combination of vector and matrix notations. Thus a 
ray will be given by its position vector: x = (zx, (¢), z2(t), 2;(t)). This is a row vector and 
a 1 X 3 matrix; the corresponding column vector will be indicated by x”. Thus the 
dot product of two vectors is simply y-z = yz’, while y’z is a 3 X 3 matrix. The unit 


3 X 3 matrix, with ones down the main diagonal, and zeroes elsewhere, will be repre- 
T 


sented by J. y° will mean y-y = yy . 
It will also be convenient to use the normalized time parameter 


u =ct. (1) 


Here c is the speed of light in a vacuum. Differentiation will be with respect to u, so 
that, for example, x’” will mean the column vector (a 3 X 1 matrix) whose components 
are dx,/du,i = 1, 2,3. Since the index of refraction, n(x, , 22 , 23), is defined as c divided 
by the speed of light in the medium, it follows that along a ray x we have 


zs? =n’. (2) 
It is at this point that the isotropic property of the medium is used. 
Of the various standard forms for the differential equations of rays, which can be 
found in Herzberger [1], Synge [3], and Luneberg [4], we will therefore use 
{n(x’?)~'/?x"}" _ (x’?)'?Tn, (3) 
where Yn means, of course, the gradient of n: 
On On on) 


ve (2, x,’ Or 





Equation (3) is usually simplified by immediate use of equation (2), but new forms 
of these equations which have direct geometric interpretations can be most readily 
derived by postponing such use. 

First note that from the chain rule for differentiation 


, 


n’ = x’-Vn. (4) 
Also, differentiating (2) yields 

x’-x’”’ = —n'*n’ = —n-*(x’- Vn). (5) 
If, then, we perform directly the differentiation of the fraction indicated in the 


left hand side of (3), and then clear (3) of fractions, a use of (2) and (4), a transposition 
and a multiplication by n yield 


n*WYn — n-'(x’-Vn)x’ = x”’ — n?(x’-x’’)x’. (6) 


*Received January 28, 1953. 
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If in (6) we transpose all matrices, we get 
{I — n?x’"x’}x"’" = {I — n*x’’x'}(n“°Vn)’. (7) 


On the other hand we could use (5) either on the left or on the right of (6), transpose 
a term and then transpose all matrices to obtain the following two forms: 


x’? = {J — 2n*x’"x'}(n* Vn)’, 
(8) 
(n~*VJn)” = {I — 2n?x’*x’}x’’’. 


It is not difficult to show, as is done in Gorn (6), that J — 2n’x’’x’ is the symmetry 
operator through the plane perpendicular to x’, while J — n’x’’x’ is the projection 
operator on to that plane. Equations (8) therefore mean that x’ and n~* Vn are sym- 
metric images through the plane perpendicular to x’, while equation (7) states that 
they have the same projection on to that plane. The three component equations of 
(7) must, consequently, be dependent, since many vectors have the same projection 
as n° Yn; (7) contains incomplete information about x” and needs, say, (2) as an 
adjunct. Either set of (8), on the other hand, contains complete information. 

At this point we can make use of the standard differential geometry of space curves 
(see e.g. Franklin [5], p. 107 on) with its concepts of curvature x, torsion 7, the natural 
coordinate vectors, namely the unit tangent T, the unit normal N, and the unit bi- 
normal B, and finally the Frenet-Serret formulae expressing the derivatives of the 
natural coordinate vectors as linear combinations of T, N, and B. Once we find the 
curvature and torsion along rays we will be able to differentiate any vector field defined 
along them by using the Frenet-Serret formulae. In particular, repeated differentiation 
of x by this means will yield power series expansions of rays. These expansions correspond 
to the one well-known in differential geometry for space curves. 

If s is the arc length measured along a ray, then s’ = n~*. Furthermore, the transpose 
of the first form of (8) is 


x’? = —2n7'(x’-Vn)x’ +n “*Vn. (9) 


From these and the curvature formula, one readily derives the well-known formulae 
for the curvature of a ray: 


«= |x’ X Vn| = —n" | Vn {sin ¥, (10) 


where y is the angle between x’ and Vn. 
The standard forms for T, N, and B, with the help of equations (5) and (9) become 


T = nx’ 
N = —nk™'(x’:- Vn)x’ + (nx) "Vn (11) 
B= ax’ X Van. 


To find a formula for the torsion of a ray, one differentiates (9), and uses the rules 
or simplifying a determinant to obtain the scalar triple product 


(x’x!’x’"") = n-*(x’ Vn [Vn]’). 


Since 


[Vn]’ = x’ Van, 
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where 


Vin = (2) (13) 
mn = \an, az;/” 


a symmetric matrix, the standard torsion formula yields 
r =x ?(x’ Vn [Vn]’) = x 7(x’ Vn x’ Vin). (14) 


Alternative forms for «x and 7 can be given in terms of the matrix 


+ 








9 me on 
OX; OX. 
an on 
y=| — 0. —-— 
O23 Ox, 
on on 
~ O22 da, . 
which is skew symmetric: Y” = —Y, and has the further properties that 


Y? = (Vn)” Vn — (Vn)’l, 
wY =wX Vn 


I 


for any row vector w, and has the characteristic polynomial \{A — (Vn)*}?. 


x? = —x’Y’x’’, 
Kr = —x’Y(V2n)x’’. 
The knowledge of « and + permits us to express the Frenet-Serret formulae for rays 
T’ = xn 'N 
N’ = —xn''T +7n7'B (15) 
B’ = —m 'N 


We now have all the apparatus needed to obtain recursion formulae for the expres- 
sion of the m’th derivative of any vector y if one knows the coefficients in y = a,T + 
boN + cB or y = dx’ + e& Vn + fo(x’ X Vn); let us begin with the first form. 


If we have 
y” =a,T + b,.N + c,B, (16) 
then differentiation and application of Frenet’s formulae (15) yield 
yo" = fal, — xn™'b,}T + {b4 + «n'a, — rnc, }N + {cl, + rn™'b,,}B. 
Therefore: 
Gos, = 2 — «n'b, , 
baer = 06 + «na, — mtn , (17) 


-1 
Cus, ™ Ce + rn b,, 
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As an example, the first equation of (11) now reads 


a4=n , 
for y = Xx, b, = 0, (18) 
c, = 0. 
Applying (17) therefore yields 
a, = —n’n 
for y=*yX, b, = xn™’, (19) 
C. = 0. 


This checks with equation (9) in view of (4) and the second equation of (11). 
Similarly, eliminating x’ from the first two equations of (11) and solving for Vn 
yields 


Mo = nn’, 
for y = Vn, bo = Mm, (20) 
Co = 0. 


Here again, then, we can apply recursion formulae (17) to give 


a, = [nn’]’ — x’, 
for y = Vn, b, = [nx]’ + xn’, (21) 
Cc; = KT. 


This formula for [Vn]’ = x’ V.n can be given an alternative form, for total deriva- 
tives of n and « can always be reduced to expressions involving partial derivatives of 
nm and x with respect to the z,; and total derivatives in x. For example, a, is the T com- 
ponent of [Vn]’ = x’ Van, hence a, = x’ V.n-T = nx’ V2nx’’. 

For the alternative form of b, , we need such an expression for x’, which, in turn, 
needs one for n’’. From the Lagrange identity and formula (10), 


ck = {/x’X Vn)? =n (Vn) —n”. (22) 

From this and the alternative expression for a, it is easy to show that 
nn’! =x! Vonx’” +n-*(Vn)? — 2n''n”. (23) 

Finally, from (23), (12), and a differentiation of (22), we find 
kK’ = —2n7'n’x — n' k(x’ Von x’”) + 07K (Tn Vn’). (24) 
Thus (21) may also be written as follows: 
a, = n(x’ V.nx’’), 

for y = Vn, b, = —nn’«'(x’ V.nx’") + 0'« (Vn Vanx’’), (21a, 


KT. 

















1953 | SAUL GORN 359 


This form for [Vn]’ permits us to return to the computation of x’ from (19) and 
(17). Naturally the same result could be obtained by using (23) and (24). In either 
case a computation yields 


: 2 ‘ T 
a; = —n*{2x? — 3n” + n(x’ V.nx’’)}, 
. -3(- ~ T - - T ~ 
for y = x, b, = —n*{5n'x + «'nn’(x’ Von x’”) — nk (Vn Van x’’)}, (25) 
C3; = Nn a Kr. 


Further application of (17) to y = x is complicated by the fact that higher order 
derivatives of n are involved beyond n”, i.e. beyond those in Vn and Y7.n; they are 
probably best handled by introducing multiple index symbols, such as are used in 
tensor calculus. At this point the matrix notation stops being convenient. So much for 
the recursion formulae of type (16). Now let us consider those of type 


y” =d,x’ te, Vn + f(x’ X Vn). (26) 
To avoid confusion let us designate the d, , e, , and f, when y = Vn by p.,q,1: 3 
if we apply (11) to (21a) we obtain 
[Vn]’ = px’ + qa Vn +7,(x’ X Vn), 


where 


Di = «°{(VYn)?(x’ Von x’”) — n'(Vn Van x’’)}, 


(27 
g, = x ?{—n'(x’ Vonx’”) + 077(Vn Vin x’’)}, 
n=T 

Differentiating (26) and using (9) with (27) yields 

Ami = di, — Qn™'n'd, + Diem + 2'T1S ms 
Cnss =e +n“d, + Gem —2 ifm; (28) 
[a + ren + {qi — 2n'n’}f, . 

As examples, if y = x in (26) and m = 1, then obviously d, = 1, e, = 0, f; = 0; 
hence, from (28), d, = —2n™'n’, e, = n-*, fz = 0, which is the same as (9). Another 
application of (28) with m = 2 yields 

d; = —2(n™'n’)’ + 4n-’n” + np, , 
for y =X, e,; = (n-*)’ — 2n-‘n’ +n “a, (29) 


f= n*r, . 
Again, if y = Vn, then d, = 0, e = 1, fo = 0, hence d, = p,,¢: = g: , f: = 7: a8 in (27). 
To obtain power series expansions of the rays 
, ] 7,2 1 (m)_m 
X = Xo + ou + 5 Xo'u + > ae UN" tire, (30) 


one can either apply (16) to y = x at u = 0 (the first four terms come from (18), (19), 
and (25) if one desires to use T, N, and B as coordinate vectors), or one can apply (26) 
and use the recursion formulae (28), as just indicated. 
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ON SUPERSONIC FLOW PAST AN OSCILLATING WEDGE* 
By MILTON D. VAN DYKE (Ames Aeronautical Laboratory, Moffett Field, Calif.) 


1. Introduction. In order to study the limitations of the linearized theory of oscillat- 
ing airfoils, Carrier’ has analyzed supersonic flow past a thick wedge which oscillates 
slightly about its apex. Considerably greater insight into the nonlinear effects of thick- 
ness is gained by generalizing the analysis to include other locations of the pivot. 

2. Analysis. Consider a wedge fixed in a supersonic stream at a Mach number above 
that at which the shock wave detaches from the apex. If now the wedge executes slight 
oscillations, the shock will remain attached, and the flow field can be found by super- 
posing a linearized acoustic field upon the nonlinear steady flow. For harmonic oscil- 
lations of frequency a,c, Carrier has shown ‘” that the velocity components and pressure 
downstream of the shock are given by 


u= IF se a.(y, + E,), 9 = aly, deans E,), (1) 
p = p.[l — vv. + Me¢.)], (2) 
where 
og = ef 2/8) $* (a, cosh v6 + b, sinh v6)J,(kr), (3) 
E=e* t-—z/M=—)y/ MB?) pm c,d (key), (4) 


and the deflection of the shock wave from its steady position is 


v= (1 — pi/pn) Te YP dS (key) (5) 


Here, U, , a2 , M, pz and p, are the flow speed, sound speed, Mach number, pressure 
and density downstream of the shock in steady flow (and p, the density upstream), 
¥ is the adiabatic exponent, ¢ the time multiplied by a, , 6’ = M’ — 1, k = c/é”, tanh 6 = 
By/z, 7? = x? — By? and & = 1 + d? — M”’. The coordinate system and X are defined 


in Fig. 1. 
*Received Feb. 18, 1953. 
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The pertinent .jump conditions across the shock wave, which can be imposed at 
the steady shock position, were shown to be equivalent to 


oT E, = (I — p:/ po) [(nay + Aae) Wy /(1 + d’) + ay,/(1 + ’)'?], (6) 
vg, — E, = (1 — pi/p2)[(az — Ane) p,/(1 +”) — Aap. /(1 + d*)"”], (7) 


Pr + Meg, = -—(1l- pi/poas[np,/(1 + y" + vi], (8) 





Fig. 1. Coordinate system. 


where a, = (2 + 4m”)/3(1 — m’), ag = (5 + m’)/6m and a, = —m(5 + m’)/3(1 — m’) 
for y = 7/5. Here m and n are the components of M normal and tangential to the un- 
disturbed shock. 

Suppose now that the wedge is pivoted at a point a distance b downstream of its 
apex, and oscillates with small angle of attack a = Re(e**‘) in some unit system. The 
boundary condition of tangent flow at the surface is* 


(gy, — E.)y-0 = —[M + ic(x — b cos e)Je", (9) 


where ¢ is the semi-vertex angle of the wedge. Following Carrier, using the generating 
function for the Bessel coefficients we find that for vy > 1 the b, are given by 


b, = tvlr? + (—1)""1/Bk + b cos er” — (—7)"'], (10) 


where 7 = 7(M + 8). 

For the special case of rotation about the apex, the summations in Eqs. (3) to (5) 
begin with vy = 1. For other pivot positions, however, a term vy = 0 must be added to 
Eq. (5) to account for the fact that as the apex of the wedge oscillates about the origin 
the attached shock wave moves with it. The question of whether corresponding modifi- 
cations are required in the series for g and F can be answered as follows. Consider the 
pivot point to be moved indefinitely far downstream, and the amplitude of angular 
oscillation correspondingly reduced, so that ultimately the wedge simply executes a 
small vertical translational (‘‘plunging’’) oscillation. Furthermore, let the frequency 
of oscillation tend to zero. In the limit, the wedge stands fixed and slightly above and 
parallel with its original position. It is clear that in this steady flow the shock wave is 
displaced from its original position, but the velocity perturbations associated with ¢ 
and £ are zero. Hence a vy = 0 mode must be added only to Eq. (5). 


’3Here an error in the tangency condition in the reference of footnote 1 has been corrected. Further- 
more, numerous typographical errors in the subsequent equations of that paper have been rectified. 
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The additional coefficient thus introduced is determined by the condition that the 

shock wave must at all times meet the apex of the wedge, which gives 
do = —b(1 — p;/p2)(1 + *)7'7(A cos € — sin ©). (11) 


Then the shock wave conditions yield three recurrence relations for determining suc- 
cessively the coefficients a, , c, and d,(v > 1) in terms of d, and the b, . In matrix form, 
with tanh 6 = 8/A, these are 


cosh v6 £ —nV1 — m? (ma, + a:)/M | | @y+1 
8 sinh v4, 0 V1 — m? (n'a, — ma) ul biel 
{—M cosh vA 0 —n 7/1 — m Qs | id, Mf 
' — cosh v6, —£ nV1— mm? (ma, + a.)/M ] ae 
+ | B sinhvA 0 ohn or (n°a, — Mae) ~ “See | 
M cosh v4, 0 nV1 — ma; id,_, but 2d, | 
—2iM cosh v4, —2ix/M 2i[m(1 — m’*)a, + n'a.) /M a, | 
+ 0 218°/M —2in[(1 — m*)a, — ma,]/M J] 
27 cosh v4, 0 2i(1 — m°)a; J d,| 
sinh v4, —sinh v4, —2iM sinh vo (byes 
- 8 cosh v6 8 cosh v6, 0 | 5,—1 == 0, (12) 
| 
—M sinh v4 M sinh v6, 27 sinh 160} lo, J 
for »y = 0,1 --- , where it is understood that the a, , b, and c, vanish for vy < 1, and the 


d, for vy < 0. Note that for vy = 1, d,_, is to be replaced by 2d, . 

3. Example: Slow oscillations. The solution can be readily converted into an 
expansion in powers of frequency. Then retaining only linear terms in. frequency shows 
that for slow oscillations the surface pressure coefficient, referred to conditions upstream 


of the shock, is 

C, = (p — p:)/40,Ui = C, — (2/Mi + yC,) [Aa + (Bb cos € + Cz)a’ /a,}. (13) 
Here C, is the value for steady flow, M, the free-stream Mach number, and @ and a 
the instantaneous angle of attack and its (true) time derivative. The coefficients A, B, C 
are 


A = —M?’na;/p B = M(n — mtan 6)a3/p, 


a 2 92 2 2 1 / et hyg2 
C=M 2u — 2M"dna;/B + (M+ N Oe mn jos/u 4. M(M + Ines (14) 
Au — Bnag B'u 











1953] JOHN W. MILES 363 


where u = n'a; — ma, . For a wedge airfoil with flat base (on which the unknown base 
pressure is assumed to be constant, though possibly time dependent) the normal force 
and pitching moment coefficients are 


Cy = (2/Mj + yC,)[2Aa + (2Bo + Cha’ sec €/az], (15) 
Cy = (2/Mi + yC,) see’ e[A(2e — lla + {Bo(2o — 1) + Clo — 2)}a’ sec e/a], (16) 
where ¢ = c cos’ ¢€/c. 

The term in Eq. (16) proportional to a represents the aerodynamic damping moment 
for slow oscillations, which tends to stabilize if it is negative. The boundary of neutral 


stability as it depends upon pivot position and free-stream Mach number is shown in 
Fig. 2 for a wedge airfoil of 5° semi-vertex angle. Also shown are the corresponding 
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Fic. 2. Boundary of neutral stability for slowly oscillating wedge airfoil of 5° semi-vertex angle. 


results from linear and second-order theory‘ (with which the present theory agrees 
when expanded in powers of ¢), which are applicable to any airfoil shape. 


A NOTE ON SUBSONIC EDGES IN UNSTEADY SUPERSONIC FLOW* 
By JOHN W. MILES (University of California, Los Angeles) 


Summary. The pressure distribution due to the unsteady motion of a wing having a 
supersonic leading edge and a subsonic trailing edge is determined by applying a Lorentz 
transformation to the corresponding result for a rectangular wing. This result, valid as 


‘Milton D. Van Dyke, On second-order supersonic flow past a slowly oscillating airfoil, J. Aero. Sci. 
20, 61 (1953). 


*Received March 2, 1953. 
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it stands for a subsonic leading edge, is then modified by the addition of an eigensolution 
in order to satisfy the Kutta condition. 

1. Introduction. It was noted in an earlier paper [1] that the general solution for a 
rectangular wing in unsteady supersonic flow could be extended to a wing having a 
straight subsonic leading edge with the aid of the Lorentz transformation. It is the 
purpose of the present note to extend this result to a subsonic trailing edge. (The corre- 
sponding steady flow problem has been treated by Lagerstrom [2] and others.) 

2. The potential equation. The notation of the present paper is in agreement with 
that of ref. 1, to which we refer by the prefix I ( ), where ( ) denotes the equation 
number therein. In the interests of brevity, the following analysis is written as a con- 
tinuation of I and repeats equations and results stated therein only insofar as appears 
necessary to define the problem at hand. 

The original coordinates (dimensionless, being referred to a characteristic length 1) 
are the orthogonal Cartesian set (x, y, z), where the free stream velocity U, is directed 
along +2. In these coordinates, the projection on z = 0 of the wing to be considered is 


bounded by 


x + mBy > 0 (2.1a) 

mBx +y = 0 (2.1b) 

—-l<mz0O0 (2.2) 

B = (M’? — 1)” (2.3) 

where x + m$y = O defines a straight, supersonic leading edge and m6@x + y = Oa 


straight, subsonic trailing edge. The definition of the leading edge is dictated by con- 
venience (so x’ = 0 there; v.7.) and may be varied at will, insofar as it remains supersonic. 
We next introduce the coordinates (£, r) by cascading Gallilean [I (2.3)] and modified 


Lorentz [I(2.4)] trarisformations as follows: 
‘ l M ] —M x 
= g' (2.4) 
T M l 0 ] at/l 


where ¢ and a represent true time and the sonic velocity, respectively. Finally, we 
introduce the oblique coordinates (x’, y’) via the true Lorentz [I(8.1)] transformation 


ctr 


zx’ | m 
y’ m ] y 
In these coordinates the leading and trailing edges are specified by x’ = 0 and y’ = 0, 
respectively. Then, taking (z’, y’, z, r) as our working coordinates, the equation for the 


velocity potential, obtained by applying the Lorentz transformation to I(2.6), reads 


Pz'z’ — dy " id Pz: _ Prt (2.6) 
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The corresponding equation for the pressure perturbation, obtained by applying 
2.5) to 1(2.12), may be written 


p = —poap'x (2.7) 


x=x(2’,y’,2,7n7 =(1- m’)~'?M(¢9,: + md,) + ¢, (2.8) 
The integration of (2.8), subject to the boundary condition ¢ = 0 at the leading 
edge, yields, 


»z’ 


d 1 — m’)'?M" | xlu, y’ — mx’ — pw), tr — (1 — m’)'?M "(2 — p)] du (2.9) 


3. The boundary conditions. The wing platform (S), as it appears in the (&, y) and 
(x’, y’) coordinates, is shown in Fig. 1. The wake aft of the trailing edge is designated 











Fic. 1. The regions S, W and R, the various coordinate axes, and the Mach lines y = +#£. 


as W and the remainder of the plane z = 0 as R. The boundary conditions in S and R are 
as in I(2.7) and I(2.10), respectively, while the boundary condition in W is dictated by 
the requirement of continuity of pressure. Thus, we obtain 


o, = —w’'(z’, y’, 7) on S(z7’ > 0, y’ > 0,2 = O+) (3.1) 
x=0 on W(x’ > 0, mz’ S y’ S$ 0,2 = 0) (3.2) 
¢=0 on R(x’ > 0, y’ < mz’, z = 0) (3.3) 


where w’, the prescribed downwash, is now exhibited as a function of the oblique co- 
ordinates. 

In the case of a non-trailing edge the region W would not exist, and (3.3) would 
hold for y’ < 0. The corresponding solution for ¢, designated as ¢“” in the following, 
and its derivatives ¢,- , and ¢, then would vanish as y’*’”” at the edge y’ = 0, while ¢,- , 
would behave as y’'”’ there. Accordingly, for m > 0 x would exhibit the expected 
singularity associated with a subsonic leading edge (cf. ref. 2). On the other hand, for 
m < 0 an eigensolution, ¢’, can be determined to cancel this singularity, thereby 
satisfying the Kutta condition (x = 0 at y’ = 0). Accordingly, we write 


¢=9o97 +o (3.4) 


and formulate two, subsidiary boundary value problems. 
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° . . . (1) ‘26 °,¢ 
As indicated in the foregoing, ¢’’ may be assumed to satisfy the boundary conditions 


¢. = —w’ on S(r’ > 0, y’ > 0,z = 0+) (3.5) 
¢” =0 on W+R(e’>0,y' S$ 0,2 = 0) (3.6) 


It follows that ¢{”’ and, therefore, x”? must vanish in S in order to satisfy (3.1). More- 
over, since ¢°”’ must vanish in R, so also must x’. Hence, we write 


x,’ =0 on S (3.7) 
x‘? =O on W+R (3.8) 


To complete the formulation, we remark that x also satisfies the differential equation 
(2.6). 

As will be shown below, the explicit determination of ¢‘’ is not required insofar as 
“) may be determined from (2.9). 


x alone is required. However, given x’, 

4. Solution. The boundary value problem of determining ¢" 
is now identical with the rectangular wing problem solved in I. A convenient form of 
the solution for our present purpose is obtained by posing the harmonic time dependence 
exp (ixr), substituting x’ and y’ for é and y in I(4.5) and introducing the trigonometric 


" , 


change of variable » = y’ + (x’ —u) cos @ cos \, whence we have 


when w’ is prescribed 


o"(x’, y’, O+, 7) = | Jolx(x’ — p)]w’(u, y’, 7) du 
I —_ ? 
—-]| du]  Jolx(x’ — u) sin 6] do (4.1) 
wr. ' 
a a2ain z 7 cos8 ° = ; é 
ary | w'[u, y’ + (x’ — uw) cos 6 cos X, rT] dA 


where the are-sine is to be replaced by +/2 when its argument exceeds unity. In the 
following, we refer to this solution in the symbolic form 


o'(x’, y’, O+, 7) = o' {w'(z’, y’, 7} 1.2 
‘ { 9 < 


To obtain x°”’, we remark that, since any derivative of ¢ also satisfies (2.6), ¢,:' {w’} 
represents a solution satisfying the boundary conditions (3.5) and (3.6) provided that 
w’ is replaced by w/. therein. Moreover, ¢‘"’{w/.} satisfies the same differential equation 
and boundary condition as ¢;'’ but differs therefrom in vanishing at y’ = 0, whereas 
¢:)’ is singular as y’'” in this neighborhood. It follows that ¢)'’{w’} — ¢°° {w}-} 
represents a solution to the differential equation (2.6) satisfying the homogeneous 
boundary conditions (3.7) and (3.8) and exhibiting the desired singularity at y’ = 0. 
Accordingly, the required solution for x” that cancels the singularity in x“ at y’ = 0 
is given by 


= —m(1 — m)'?M[g<? {w’} — o'? {wi 4] (4.3) 
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Substituting from (4.1), we have 


x’ cos~'(y’/(2’—z)] 
x? m(1 — m’*)"'? My’ “_ [ du [ J o[{x(a’ — pw) sin 6)] 


“0 0 
(4.4) 
9 ; 
= {{(x’ — p) cos 06 — y’|'*w’[u, (x’ — pw) cos 0 — y’, r]} dé 
€ 
The final solution for the pressure distribution on the upper surface of the wing, 
obtained by substituting 6‘” in (2.8) and adding x from (4.3), is given by 


x(x’, y’, OF, 7) = (1 — mm?) OMe! {[w’} + mo™ {wt.}] + 62? {w’} (4.5) 


We remark that the results (4.3) and (4.5) are not restricted to harmonic time dependence, 
since they are valid for all frequencies. 

In the case of a subsonic leading edge (m > 0), it would be necessary only to replace 
¢'’ {wi.} by dy’ {w’} in (4.3). 


REFERENCES 
1. J. W. Miles, A general solution for the rectangular airfoil in supersonic flow, Q. App]. Math., 11, 1-8 
(1952). 
2. P. A. Lagerstrom, Linearized supersonic theory of conical wings, J. P. L. Prog. Rep. 4-36 Pasadena, 
(1947); reprinted as NACA TN 1685 (1948), 90-98. 


NOTE ON THE MEAN SQUARE VALUE OF INTEGRALS IN THE 
STATISTICAL THEORY OF TURBULENCE* 


By C. C. LIN (Massachusetts Institute of Technology) 


1. In the statistical theory of homogeneous isotropic turbulence, it is sometimes 
of interest to evaluate the mean square value of certain integrals, such as the pressure 
fluctuation over a sphere. The purpose of the present note is to give such an evaluation 
for integrals over a sphere and for similar integrals over spaces of other dimensions. 
The analysis shows that the final answer can be interpreted in terms of dimensional 
arguments; provided the length scale used is the geometrical mean of the scale of turbu- 
lence and the linear scale of the region over which the integral is taken. The results could 
be applied to the problem of the noise generated by turbulence. 

2. Consider, for definiteness, the pressure fluctuation over the surface of a sphere. 
Extension to the study of other quantities can be easily made. Let the integral be 
denoted by 


[= | vas, (1) 


where p is the pressure fluctuation at a point P, and the surface integral is extended over 
a sphere of radius a. We may also write 


f= / p’ ds’, (2) 





*Received April 20, 1953. 
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where p’ is the pressure fluctuation at a point P’, and the integral is extended over the 
same surface of the sphere. Thus, if we calculate J° by multiplying (1) and (2), we obtain 


(I*) = [ as [ ip’) dS’, (3) 


where ({ ) enclosing a quantity denotes a statistical average. 
In isotropic turbulence, the statistical correlation (pp’) depends only on the relative 
position of the two points P and P’. Thus, the integral 


J = |’ dS’ (4) 


is independent of the position of the point P, and we have 
(1°) = 4na’J. (5) 
To evaluate the integral J, we may take the point P as the origin of a system of 


spherical coordinates. The element of area can be represented very simply if we con- 
sider zonal surfaces at a distance r from the origin. In fact, 

dS = 2zxr dr. (6) 
Thus, 


J = 2r(p’) | P(r)r dr, (7) 


“0 


where (p’) is the mean square value of the pressure fluctuation, and @(r) is the correla- 
tion coefficient for pressure. 
For very small spheres, we have the approximation 


P(r) ~ 1 for 0<r< 2a. (8) 


Then J may be approximated by 


Jo = 4na°(p’) (9) 
and (7’) may be approximated by 
(I*), = (4a’)*(p’) (10) 
This also follows, if we write 
I = 4ra‘p (11) 


for very small spheres. 
For large spheres, J may be approximated by 


J. = 2p’) | @(r)r dr. (12) 


#0 


This may be written as 


Jo = 4x(p yl, (13) 
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where / is a scale of turbulence defined by 


2° = [ @(r)r dr. (14) 
Then 
(I*) = (4na*)(4nT°)(p’). (15) 
The ratio (J*)/{I*), is equal to 
.. (4) 
(P), ~ \a (16) 


In the limit a +, this approaches zero. This is to be expected since [/4za’, being the 
average value of J over the sphere, approaches zero when the sphere increases indefi- 
nitely in size. 

In many cases, we are interested in averages over spheres of a size larger but com- 
parable to the scale of turbulence; then (15) and (16) may be used as suitable approxi- 
mations. 

3. From the above arguments, it is obvious that the order of magnitude of I’) is 
independent of the detailed shape of the surface under consideration. Thus, we may 
write 

(I?) = kp’\PL’, (17) 
where k, is a constant, depending somewhat on the shape, / is a scale of turbulence, and 
L is a typical linear scale of the surface. 

The formula (17) may be easily generalized to other dimensions. In general, we have 

(2) = k,(p’)(1L)", (18) 


for an integral 
I, = [ dr, (19) 


over an n-dimensional space. Thus, (18) can be obtained from (19) by dimensional 
arguments provided the volume integration is associated with the nth power of the 
length scale (lL)', i.e., the geometrical mean of the scale of turbulence and the linear 
scale of the region over which the integral is taken. 

4. Some care must be exercised in estimating such integrals, if the divergence theorem 
can be used to convert them to integrals of lower dimensions. Consider, for instance, the 
integral 


ou 
= ea dr (20) 


K 
over a sphere. If we apply (18) directly to (20), we obtain 
(K*) = k,((du/dx)*)(IL)*”?, 
or 
(K?) ~ 2.) (21) 
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where \ is Taylor’s micro-scale of turbulence. However, this is actually incorrect. If 
we write (20) in the form 


K = [ un, dS, (22) 


then it becomes evident that 
(K*) ~ (u’\(IL)’. (23) 


The ratio between the two estimates (21) and (22) is IL/d*. It is clear that the direct 
application of the volume formula over-estimates the integral because of the small scale 
introduced by the differentiation process. Such occurrences are frequent in the statistical 
theory of turbulence. For example, the divergence of the Reynolds stress 7;; = —p(u,;u;) is 


Or re] ‘ , Ou; 
. = — };—p(u,u;)} = —p\u; = /. 
OX Oz; ' p “nes Ox; 


| 
Based on the left-side, the estimate is 


OT; T 
—— fay 


Ox; l 


while the right-side would give the (incorrect) estimate 


Or; T 


Ox; 


In all such cases, the lower estimates are to be taken. 


BOOK REVIEWS 


Gasdynamik. By K. Oswatitsch. Springer-Verlag, Wien, 1952. viii + 456 pp. $18.60. 


To write a comprehensive treatise on a rapidly expanding field of knowledge is not an attractive 
task: there is always the possibility that new developments may soon make the treatment appear dated 
or even incomplete. On the other hand, the lack of such a treatise may seriously impede the recruiting 
of new scientific talent, because newcomers will find it increasingly difficult to work their way through 
numerous important papers written in various languages and scattered over a great number of technical 
periodicals. In producing this comprehensive treatise on the dynamics of compressible fluids, the author 
has therefore rendered a significant service to all interested in the development of this branch of mechanics 
of continua. 

Chapter I contains the necessary thermodynamic background. Chapters II and III are concerned 
with steady and unsteady flows in one dimension. The fundamental integral theorems are established 
in Chapter IV. These integral theorems remain valid in the presence of shocks; the differential equations 
of motion are readily derived from the integral theorems. Mechanical similarity is discussed, and various 
vortex theorems are presented. Chapter V illustrates the application of the integral theorems to technical 
problems. Chapter VI is devoted to the general equations for steady inviscid flow and to exact particular 
solutions of these equations (Prandtl-Meyer flow, axially symmetric conical flow, transformations of 
Molenbrock and Chaplygin, linearization of Prandtl and Glauert). Chapter VII is concerned with 
steady subsonic flows (plane or axially symmetric). In particular, the methods of Krahn, Janzen-Ray- 
leigh, Karman-Tsien, and Ringleb are discussed. Chapter VIII is devoted to steady supersonic flows 
in two dimensions (slightly disturbed parallel flow, shocks and their interaction, method of character- 
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istics). Chapter IX surveys the available results regarding transonic flow. Chapter X is concerned with 
particular steady and unsteady three-dimensional flows (lifting surfaces, conical flows, delta-wing, 
decelerated wedge). The influence of viscosity is briefly discussed in Chapter XI, and a survey of ex- 
perimental techniques is given in Chapter XII. 

W. PRAGER 


An introduction to linear programming. By A. Charnes, W. W. Cooper, and A. Henderson. 
John Wiley & Sons, Inc., New York, and Chapman «& Hall, Ltd., London. ix + 74 
pp. $2.50. 


The first part of the volume contains an economic introduction to the theory of linear programming. 
The fundamental] concepts are presented in connection with a specific problem (“nut mix problem’’), 
and the mathematical argument is kept as simple as possible. The second part is devoted to the mathe- 
matical theory of linear programming. 

While all concepts with which the non-mathematical reader is not likely to be familiar are fully 
explained, the second part, nevertheless, requires greater mathematical experience from the reader 
than the first part. Obviously, there is some duplication of the mathematical work in the two parts, 
but this is entirely justified in a work which is aimed at such a diversified group of readers. 

W. PRaAGER 


A selection of tables for use in calculations of compressible airflow. Prepared on behalf 
of the Aeronautical Research Council by the Compressible Flow Tables Panel. 
The Clarendon Press, Oxford, 1952. viii + 143 pp. $8.00. 


The panel consisted of L. Rosenhead (Chairman), W. G. Bickley, C. W. Jones, L. F. Nicholson, 
C. kK. Thornhill and R. C. Tomlinson. 

This book of tables is the first of two volumes. A complementary volume consisting of graphs will 
be published in the near future. Among other tables, there is a series of tables labeled isentropic flow 
tables, applying to the steady isentropic flow of air; a series of tables, labeled characteristic tables, 
applying to steady isentropic supersonic flow of air in two dimensions or in three dimensions with axial 
symmetry; and a series of tables, labeled shock tables, applying to normal and oblique shocks in air. 

P. CHIARULLI 


Advanced mechanics of materials. By F. B. Seely and J. O. Smith. John Wiley and Sons, 
Inc., New York, and Chapman and Hall, Ltd., London, 1952. Second edition. xvii + 
680 pp. $8.50. 


This second edition of the book first published in 1932, and then by Seely only, is essentially a new 
book. It is written for advanced undergraduate and first year graduate students in engineering. Special 
emphasis is laid on the engineering evaluation of analyses of problems. Results are sometimes quoted 
from cited references, e.g., in the discussion of thin plates, and thereby the mathematics is kept at an 
elementary level. The discussion of energy methods and of plastic instability of struts is commendably 
clear. Problems and references are provided. 

The book is in six parts [Preliminary considerations. Special topics on the strength and stiffness of 
members subjected to static loads. Localized stress—stress concentration. Energy methods. Influence 
of smal] inelastic strains on the load-carrying capacity of members. Introduction to instability—buckling 
loads. ] with three Appendixes [A brief introduction to the mathematical theory of elasticity. The elastic 
membrane (soap film) analogy for torsion. Properties of an area. ] 


H. G. Hopkins 
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Theory of elasticity and plasticity. By H. M. Westergaard. Harvard University Press, 
Cambridge Mass., and John Wiley & Sons, Inc., New York, 1952. xii + 176 pp. $5.00. 


Professor Westergaard’s original plan was to write a definitive work on two and three-dimensional 
elasticity and to include plasticity. A chapter on history covers analytical and experimental develop- 
ments in elasticity, photoelasticity, plasticity, analogies, structures, and mechanics of materials. Realiz- 
ing he would not have time to finish, Dr. Westergaard made a great effort to complete the part which 
has been published posthumously. An interesting feature is the emphasis placed upon the displacements, 
strains, and the strain functions from which they may be derived rather than upon stress. The approach 
throughout employs vector notations and descriptions wherever possible. Galerkin vector solutions of 
the problems of Kelvin, Boussinesq, Cerruti, and Mindlin are discussed in detail as is the author’s own 
method of the twinned gradient 

D. C. Drucker 


Mechanics of vibration. By H. M. Hansen and P. F. Chenea. John Wiley & Sons, Inc., 
New York and Chapman & Hall, Ltd., London, 1952, xiii + 417 pp. $8.00. 


This book is intended as a textbook for both undergraduate and graduate engineering students. 
In the authors’ words, it is to be a textbook that covers the all-important basic principles in a 
ble for a student who has had nothing more than an elementary course 


thorough fashion and vet is suita 
standard instruction in mathematics offered to engineering students today’’. This 
reviewer is of the opinion that the authors have accomplished their purpose in large measure insofar as 
idents are concerned, but that the book is not well suited as a text for graduate students. 


in dynamics and the 


undergraduate sti 
Part I treats at great length a) free vibrations without damping, b) harmonic foreed vibrations 


without damping, the steady state solution of harmonic forced vibrations with damping. of linear 


systems having one degree of freedom. 
Part II treats systems of several degrees of freedom. In chapter five there is a brief discussion of 
ions and of the solution of the equations of small oscillations. In the reviewer’s opinion 


Lagrange’s equat 
the book because the all-important subject of the fundamental theory 


this is the weakest chapter In 
of small oscillations is discussed only very sketchily. Evidently the authors were prevented by space 
limitations from giving a fuller account of this theory. This reviewer believes that a condensation or 
omission of other topics, for the sake of a more thorough presentation of this theory, would have been 
desirable, particularly since the book is intended to be suitable for graduate as well as undergraduate 
students. This chapter contains a few inaccuracies 

Chapter six is the most valuable portion of the book. It contains an exhaustive discussion of the 
mobility method for the analysis of certain types of systems of several degrees of freedom. The treat- 
ment of various lever systems, branched systems and multimass systems will be very useful to anyone 
faced with the investigation of problems to which this method applies. Since the method is based on 
systems, 


mpedance which is most conveniently developed in connection with electrica 


the concept o j 
the reviewer believes that it is unfortunate that no electrical systems are mentioned. 

Part II closes with a discussion of two methods of approximate calculation of the frequencies, 
viz., Holzer’s and Graeffe’s method. 

Part III, entitled “Special Topics’’, consists of three chapters, one on continuous systems, one on 
vibrations of transient character, and one on non-linear vibrations. These are well written, but, being 
relatively short, barely touch on the subjects involved. 

There is a large and remarkable collection of useful problems with answers. 

The book as a whole is clearly written. The mathematical tools used are all elementary. This re- 
viewer was surprised at the almost complete lack of literature references. 


G. W. Morcan 
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Elasticity in engineering. By Ernest E. Sechler. John Wiley & Sons, Inc., New York, 
and Chapman & Hall, Ltd., London, 1952. ix + 419 pp. $8.50. 


This book is intended to bridge the gap between texts on strength of materials and engineering 
treatises on the theory of elasticity. It consists of three parts: Fundamental Equations and Analysis 
Methods, Engineering Problems in Stable Structures, Engineering Problems in Instability. The subject 
matter was evidently selected with a view toward the interest of structural engineers, in general, and 
aeronautical engineers, in particular. 

The emphasis on theory and understanding is rather light. A few examples may suffice to illustrate 
the point. The local transformation theory for a state of stress is treated on the assumption of a homo- 
geneous stress field and under ‘‘neglect’’ of the body forces. In a chapter on stress functions Maxwell’s 
stress functions are unearthed and applied exclusively to the representation of a uniform stress field. 


Under the heading “Uniqueness of Energy Solution” the author claims to prove that for an elastic 
body in equilibrium under a single concentrated force the displacement at the point of application of 
the force is unique. This reviewer is unable to follow the argument. 


E. STERNBERG 


Mechanics. By A. Sommerfeld. Translated from the German by M. O. Stern. Academic 
Press, Inc., New York, 1952. xiv + 289 pp. $6.50. 


This book constitutes volume 1 of the author’s Lectures on theoretical physics. It is concerned with the 
mechanics of systems of a finite number of degrees of freedom. (Systems with an infinite number of 
degrees of freedom are treated in vol. 2 entitled Mechanics of deformable bodies an English translation of 
which was published in 1950.) The following chapter headings will give an idea of the scope of the book: 
Mechanics of a particle—Mechanies of systems, principle of virtual work, and d’Alembert’s principle— 
Oscillation problems—The rigid body—Relative motion—lIntegral variational principles of mechanics 
and Lagrange’s equations for generalized coordinates—Differential variational principles of mechanics— 
The theory of Hamilton—Problems. There is no need to comment here on the merits of Sommerfeld’s 
exposition; the translation is well done. 


W. PRAGER 


Grenzschicht-Theorte. By H. Schlichting. G. Braun, Karlsruhe, 1951. xv + 481 pp. 
$10.00. 


The key papers on boundary layer theory are scattered over many technical journals some of which 
are not readily accessible. Moreover, development has been very intense in this field for the last three 
decades and is still continuing at a rapid pace. These circumstances make a comprehensive treatment 
such as Professor Schlichting’s book particularly valuable, for the research worker as well as the new- 
comer to the field. The subject matter is organized into four parts: Basic laws of viscous flow—Laminar 
boundary layers—Transition from laminar to turbulent flow—Turbulent boundary layers. In addition 
to the material indicated by the heading, the first part contains a valuable review of exact solutions of 
the Navier-Stokes equations. The second part presents the theory of the laminar boundary layer for 
plane steady flow, for axially symmetric, and three-dimensional flow, and for unsteady, in particular 
oscillatory flow. Temperature boundary layers and boundary layers in compressible flow are also dis- 
cussed. The third part is rather brief, covering only thirty-one pages. The discussion of turbulent flows in 
the fourth part is based on the older hypotheses regarding mixing length or similarity rather than the 
more recent statistical concepts. The exposition is clear throughout and easy to follow. 

W. PRaGER 
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Advanced dynamics. By E. Howard Smart. Macmillan and Co., Ltd., London, 1951. 
Vol. I, xi + 419 pp. Vol. II, xi + 420 pp. $12.00 per set. 


One outstanding feature of this two-volume work is its tremendous scope. No previous knowledge 
of mechanics is assumed, and the subject of the dynamics of a particle and of rigid bodies is developed 
slowly, methodically, and in great detail from the very beginning through to the more advanced topics 
such as the Hamilton theory. Many topics, which in many books on mechanics are given only a cursory 
treatment, appear here in great detail. The text abounds in examples worked out at great length. There 
are extensive exercises for the reader, many of the problems being quite interesting. As for prerequisites, 
a first course in differential equations should, for the most part, suffice for volume I and the first half 
of volume II. However, for the last half of volume II considerable more mathematical maturity seems 
advisable. 

Volume I deals with the dynamics of a particle, almost all of the material being on two-dimensional 
topics. The kinematics of a particle is first introduced, then Newton’s laws of motion. There follow 
extensive treatments of the rectilinear motion of a particle, as well as the plane motion of a particle 
investigated by the use of acceleration components in the directions of rectangular cartesian coordinates 
and plane polar coordinates. Included here are broad treatments of the ballistic problem and central 
orbits. Among other topics are the motion of a particle on a curve or surface and the motion of chains. 
Volume I concludes with a chapter on the motion of a particle in three dimensions. 

Volume II deals with the dynamics of a rigid body, and with certain advanced topics. Two-di- 
mensional motions are first considered at some length. Included here is an extensive treatment of im- 
pulses. There is one chapter on three-dimensional motions. About the last half of volume II is occupied 
with more advanced topics including the top and gyroscope, Lagrange’s equations of motion, the Hamil- 
ton theory, and the theory of small oscillations of conservative systems. 

The reviewer has but two criticisms to make, both of which may be matters of opinion. First, the 
title of the work seems a bit misleading, since the subject of dynamics is developed in this work right 
from the beginning, and only the last half of volume II deals with what is commonly called “advanced 
dynamics’’. Secondly, almost no use of vectors is made. The notion of a vector is introduced after it 
has been thoroughly motivated by the velocity and acceleration of a particle. But the only operation 
introduced from the algebra and calculus of vectors is addition. Very many derivations and equations 
are thus expressed in terms of components, which is relatively awkward and space consuming. Also, no 
special notation is used for vectors, so that it is sometimes difficult to ascertain whether vector or scalar 


character is implied. 
G. E. Hay 


Numerical solution of differential equations. By W. E. Milne. John Wiley & Sons, Inc., 
New York, and Chapman & Hall, Ltd., London, 1953. xi + 275 pp. $6.50. 


As is stated in the preface, the book attempts to acquaint the reader with the principal techniques 
for the numerical solution of ordinary and partial differential equations. 

Chapters 1 (Introduction), 3 (Analytical Foundations), and 9 (Linear Equations and Matrices) 
contain useful background material not directly connected with numerical techniques of integrating 
differential equations. Elementary methods of solving ordinary differential equations are discussed in 
Chapter 2 (point-slope formula, trapezoidal formula). Chapter 4 is devoted to methods based on numer- 
ical integration and Chapter 5 to the method of Runge-Kutta and methods involving higher derivatives. 
Systems of ordinary differential equations are treated in Chapter 6, and two-point boundary value 
problems in Chapter 7. Chapter 8 deals with explicit methods for parabolic and hyperbolic partial differ- 
ential equations, and Chapter 10 with implicit methods for elliptic equations. Numerical methods of 
obtaining characteristic numbers are presented in Chapter 11. Appendices are concerned with round-off 
errors, large-scale computing machines, and the Monte-Carlo method. 

The exposition is careful and clear. The detailed numerical examples will enable the practical com- 
puter to acquaint himself with the characteristic features of the various methods, Reflecting the consid- 
erable practical experience of the author, the book is a most welcome addition to the literature on this 


important field. 
W. PRAGER 
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Vorlesungen tiber die Theorie der Integralgleichungen. By I. G. Petrovskij. Translated 
from the Russian by R. Herschel. Physica-Verlag, Wiirzburg, 1953. 100 pp. $7.80. 


The little volume contains an excellent exposition of the Fredholm-Hilbert-Schmidt theory of linear 


integral equations. 
W. PRAGER 


Thermionic vacuum tubes and their applications. By W. H. Aldous and Sir Edward 
Appleton. John Wiley & Sons, Inc., and Methuen & Co. Ltd., London, 1952. vii + 
160 pp. $2.00. 


This monograph treats the basic theory and applications of high vacuum tubes including velocity 
modulated tubes (klystrons, etc.). Written in a clear and concise manner it serves the authors’ purpose 
as an introduction to electronic circuits while the chapter on limits to amplification is worth reading 


for the more advanced worker. 
8S. L. Levy 





International Congress of Mathematicians 1954 


The International Congress of Mathematicians 1954 will be held in Amsterdam from 
September 2nd to September 9th under the auspices of ,,Het Wiskundig Genootschap”’ 
(The Mathematical Society of the Netherlands). It is the sincere hope of the ,, Wiskundig 
Genootschap” that the Congress 1954, which will be open to all mathematicians from all 
parts of the world, will be a fertile international gathering. 


The Organizing Committee has invited a number of outstanding mathematicians to 
deliver one-hour addresses, hoping that in this way a survey of the recent development 
in the whole field of mathematics may be furnished. 

There will be seven sections, viz: 

Algebra and Theory of Numbers. 

Analysis. 

Geometry and Topology. 

Probability and Statistics. 

Mathematical Physics and Applied Mathematics. 
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Logie and Foundations. 
Philosophy, History, and Education. 


“I 


In each of these sections half-hour addresses will be delivered by experts on invitation 
of the Organizing Committee. Moreover short lectures will be given by members of the 
Congress who have applied beforehand to the Organizing Committee. The time allotted 
for each short lecture will be 15 minutes. It will depend on the number of these short 
lectures whether and how the sections will be divided into subsections. 


The Organizing Committee is planning several entertainments and also a number of 
interesting excursions. 

There will be two categories of membership in the Congress: 
regular members (members) who will be entitled to participate in the scientific and social 
activities of the Congress and to receive the Proceedings of the Congress, and 
associate members who, accompanying regular members of the Congress, do not participate 
in the scientific programme nor receive the Proceedings, but will be entitled to many of 
the privileges of membership. 


The fees to be paid by members and by associate members have not yet been definitely 
fixed but presumably they will not exceed the amount of fifty guilders (about $14.—) for 
members and twenty guilders ($5.50) for associate members. 

Those who wish to attend the Congress are requested to communicate their name 
(with degrees, qualifications, etc.) and full address to the Secretariat as soon as possible. 
They will receive a more detailed communication which will be sent out in the course of 
1953. 


Amsterdam, 2d Boerhaavestraat 49 The Organizing Committee 


























